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1 . INTRODUCTION 

Whenever  a simulation  problem  is  to  be  programmed  for  hybrid  conputation,  the 
first  question  to  be  answered  is  how  should  the  problem  be  partitioned  between 
the  analogue  and  digital  portions  of  the  hybrid  conputer?  One  naturally  leans 
towards  a scheme  which  is  convenient  for  programming,  and  which  allows  easy  change 
of  system  parameters  and  equations  during  a series  of  runs.  The  complicating 
factor  is  that  due  to  the  effects  of  time  delays  caused  by  calculations  on  the 
digital  computer,  and  the  sample-and>hold  operations  of  the  interfaces  between 
the  two  computers,  the  chosen  configuration  may  not  be  mathematically  stable. 

Such  time  delays  are  always  present  in  a true  hybrid  simulation,  and  it  would  be 
useful  to  know,  in  advance  of  actually  programming  a configuration,  what  time 
delays  (i.e.  digital  coaputer  step-length)  can  be  tolerated  within  the  proposed 
partitioning  scheme,  if  the  simulation  is  to  run  in  real  time.  A considerable 
amount  of  time  and  effort  may  be  wasted  in  partitioning  and  programming  a problem 
on  a trial  and  error  basis. 

Reference  1 provides  a method  of  analysis  which  will  allow  the  user  to 
deteraune,  in  advance  of  programming,  the  stability  or  otherwise  of  a given  hybrid 
siiaulation  for  various  step- lengths,  provided  the  user  can  find  a linear  approx- 
imation to  the  response  of  the  system  under  consideration.  One  of  the 
recoamiendations  of  reference  1,  viz,  that  the  stability  analysis  procedure 
addressed  in  the  report  be  implemented  as  a software  package,  has  now  been 
completed. 

The  material  in  this  report  is  a description  of  the  software  package,  including 
information  required  by  a potential  user.  Section  5 deals  with  two  simulation 
problems,  and  shows  the  method  by  which  the  linear  representations  of  the  systems 
arc  defined  in  a form  suitable  for  input  to  the  software  package.  One  of  these 
cxaiiq>les  is  further  considered  with  a much  increased  speed  of  response.  'Hiis 

sutisequently  allows  the  question  of  stability  of  the  modifod  system  to  be  settled 
by  theoretical  analysis.  This  example  provides  a useful  test  case  for  the 
software  package.  In  addition,  one  of  the  applications  of  Section  5 has  been 
checked  by  hybrid  simulation  giving  good  agreement. 

2.  BRIEF  OVERVIEW  OF  THE  HYBRID  PROBLEM 

Once  a user  of  the  hybrid  computer  has  decided  upon  a partitioning  schesie,  he 
needs  to  answer  two  further  questions.  The  first  is  how  long  will  the  digital 
portion  of  the  hybrid  coBg>uter  take  to  coiqpute  those  parts  of  the  simulation 
problem  assigned  to  it?  The  time  for  computation,  h,  provides  a minimum  value 
for  the  simulation  step-length.  The  scheme  must  be  stable  for  at  least  this 
step- length,  and  since  there  are  other  delays,  such  as  sanple-and-hold  delays 
associated  with  hybrid  computation (ref .2) , it  is  clear  that  the  step-length  will 
need  to  be  at  least 


h ■ h ♦ fih  s 
s 

where  5h  represents  the  accumulation,  per  step,  of  these  additional  delays. 

The  software  package  provides  a quick  and  flexible  means  of  determining  the 
answer  to  the  second  problem,  viz,  will  the  proposed  partitioning  scheme  be 
stable  for  the  minima  allowable  step- length  h^?  In  practice,  the  second 

question  is  the  one  to  be  addressed  first. 

Once  a system  has  been  linearized  it  is  possible,  for  little  further  effort, 
to  examine  a nuriar  of  different  partitioning  schemes  with  repeated  use  of  the 
software  package.  Only  after  careful  selection  of  a scheme  from  the  possible 
contenders  is  it  time  to  address  the  first  question  posed  above,  and  then  usually 
only  by  specialists  familiar  with  the  particular  hybrid  installation.  It  may 
turn  out  that  the  only  way  to  achieve  an  accurate  enough  conputation  time  is  to 
program  the  equations  assigned  to  the  digital  co^>uter  and  run  then  on  the  digital 
portion  of  the  hybrid  installation.  Saag>le-and-hold  delays  must  be  added  to  the 
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computation  time,  before  the  user  is  in  a position  to  establish  whether,  h^,  the 

minimum  step- length  for  the  hybrid  simulation,  is  within  the  step- length  for 
mathematical  stability. 

The  method  of  determining  the  stability  of  a system  reduces  to  determining 
the  stability,  or  otherwise,  of  integration  processes.  As  such,  the  method  is 
applicable  only  for  certain  classes  of  systems,  namely  those  which  can  be  linear- 
ized about  a current  state  of  motion.  However,  this  can  be  realized  for  a great 
many  systems  of  interest,  since  it  is  not  required  to  provide  an  accurate 
simulation  of  system  behaviour,  only  to  approximate  the  speed  of  response 
sufficiently  well  for  the  purposes  of  stability  analysis. 

^pendix  I,  and  further  analysis  given  in  reference  1,  show  that  the  stability 
requirements  for  a given  scheme  of  partitioning  can  be  derived  in  a straight- 
forward manner.  The  ease  with  which  this  can  be  done,  however,  belies  the  care 
needed  in  the  selection  of  an  acceptable  scheme.  For  although  the  stability 
analysis  procedure  described  in  Appendix  I may  show  that  a particular  configur- 
ation is  stable,  the  stability  may  be  only  marginal.  Since  we  are  employing  a 
linear  approximation  to  the  real  simulation  model,  the  solution  for  the  latter 
may  therefore  turn  out  to  be  unstable.  In  particular,  it  should  be  recognized 
that  the  partitioning  schemes  under  consideration  allow  in  principle  for  the 
part  of  the  simulation  on  the  analogue  computer  to  be  unstable,  even  though  the 
overall  simulation  is  not. 

It  was  shorn  in  reference  1 that  relatively  poor  stability  results  whenever 
integrations  are  carried  out  digitally  and  derivatives  are  evaluated  on  the 
analogue  computer,  or  derivatives  are  evaluated  digitally  and  the  remainder  of 
the  simulation  is  carried  out  on  the  analogue  computer.  The  strong  inference 
is  that  because  of  the  relatively  poor  stability  of  either  division,  one  should, 
wherever  possible,  integrate  derivatives  on  the  same  portion  of  the  hybrid 
computer  as  they  are  calculated.  In  particular,  "fast"  loops  should  be 
simulated  on  the  analogue  portion  or  if  this  is  not  possible,  entirely  on  the 
digital  computer. 

This  does  not  mean  that  one  should  rule  out  alternative  schemes.  In  fact, 
the  recommendation  is  that  the  stability  analysis  procedure  described  in 
Appendix  and  inp lamented  in  the  software  package,  be  used  to  determine  the 
stability  of  partitioning  schemes  devised  on  the  basis  of  convenience  to  the  user. 
However,  if  there  are  no  such  requirements,  the  principles  described  in  the  last 
paragraph  should  be  followed. 


3.  REQUIREMENTS  OF  THE  SOFTWARE  PACKAGE 

This  Section  defines  the  partitioning  problem  mathematically,  and  lists  the 
procedures  required  of  a software  solution  for  analysing  stability  of  a system. 

3. 1 Problem  definition 

Following  reference  1,  we  consider  the  problem  of  solving  on  the  hybrid 
computer  the  set  of  differential  equations 

4 = Ai  (1) 

which  are  obtained  by  linearization  of  the  defining  equations  of  the  system 
to  be  simulated. 

Figure  1 shows  the  general  partitioning  scheme  for  equation  (1),  in  which 
integrations  and  evaluation  of  derivatives  can  be  performed  digitally  or  on 
the  analogue  without  restriction.  The  analysis  of  the  scheme  is  provided  in 
Appendix  I,  with  further  explanation  in  reference  1. 

The  equations  can  be  rewritten  as 
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1 - A,x  + Aii  ♦ DiX,  ♦ Dji 

i = Ajjr  ♦ A4i  + DjX,  + D4i 


where 


Further,  x is  evaluated  on  the  analogue  portion  and  X is  evaluated 
digitally,  and  the  vectors  AiX*  AsX,  A4X  are  computed  on  the  analogue 

portion,  while  DiX»  Dii.  DsJC*  determined  digitally. 

In  practice,  for  all  but  the  most  complex  partitioning  schemes,  some  of 
the  matrices  Aj , A3,  A*,  Dj , Da  and  D4  will  be  null,  and  the  others, 
including  A|  and  D|  will  be  sparse. 

3.2  Four  steps  required  of  the  software  package 

The  main  steps  in  the  procedure  for  analysing  the  stability  of  the  general 
partitioning  problem  areCsee  ref.l):- 

(1)  Specify  the  matrices  Ai  , ...  A4,  Di , ...  D4 . 

(2)  Determine  a fundamental  matrix  for  Ai  (following  procedures  described 
in  Appendix  I of  reference  1)  by  two  methods,  providing  a user  option 
for  either,  namely 

(i)  the  direct  approach  using  eigenvalue  theory,  or 
(ii)  numerical  integration  over  (o,h)  for  a given  h. 

(3)  For  a given  value,  or  values  of  h,  determing  G(h)  and  ensure  that  G(h) 

corresponds  to  a stable  analogue  simulation,  and  Q(h),  the  latter  by 
numerical  integration  over  (o,h);  G(h),  Q(h)  are  described  in 

Appendix  I. 

(4)  Construct  D and  determine  its  eigenvalues,  for  each  desired  value  of 
h,  where 

G(h)  Q(h)  (A,  + D,)  I 

h(A3  ♦ D3)  4 h D4  I 

Q(h)  D,  St  I 

Si  h A4  I 


Some  of  the  procedures  in  steps  (1)  to  (4)  above  can  be  carried  out  using 
standard  library  subroutines.  Further,  the  process  can  be  repeated  to 
cover  a range  of  values  of  step- length,  h. 

3.2.1  TWO  ^iproaches  for  step  (2),  and  merits  of  each 

Initially,  the  direct  method  for  calculating  the  fundamental  matrix 
4 of  A| , was  chosen  as  being  superior,  because  it  gives  an  exact 
solution,  whereas  the  numerical  integration  technique  gives  an 
approximate  solution.  However,  it  was  necessary  to  cater  for  both 
approaches  in  the  software  package,  since  the  properties  of  the 
eigenvalues  of  A|  may  restrict  the  use  of  one  or  other  method.  A 
method  of  analysis  was  derived  and  is  described  in  Appendix  II,  which 
allows  a choice  of  step- length  for  the  numerical  integration,  based  on 
stability  and  accuracy  of  the  solution.  In  fact,  the  results  for  both 
methods  for  obtaining  <Kt)  did  not  differ,  for  the  purpose  of  determining 
stability. 
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It  appears  that  neither  method  has  any  particular  merit  over  the 
other,  except  when  the  properties  of  the  eigenvalues  of  Ai  restrict  the 
use  of  either  method. 

It  is  suggested  that  the  numerical  integration  method  be  used  as  a 
first  approach.  Regardless  of  the  choice  of  method  for  computing 
‘l>(t),  the  software  package  calculates  the  eigenvalues  of  Ai  . 

3.2.2  Further  comments  on  step  3 

It  is  possible  to  devise  a partitioning  scheme  in  which  the  analogue 
portion  is  unstable.  This  form  of  partitioning  is  considered 
undesirable  and  the  software  package  does  address  this  question.  In 
step  (3)  the  requirement  states  that  G(h),  and  its  stability,  be 
examined.  G(h)  corresponds  to  the  analogue  part  of  the  simulation. 
This  reduces  to  determining  that  the  eigenvalues  of  G should  not  be 
greater  than  unity  in  modulus.  This  will  be  true  if  and  only  if 
there  are  no  eigenvalues  of  A|  with  positive  real  part.  This 
constraint  is  incorporated  into  the  software  package  output  to  inform 
the  user.  However,  the  package  will  continue  to  determine  the 
stability  of  the  overall  simulation  using  the  direct  method  to  compute 
♦(h) , but  the  user  is  reminded  that  such  stability,  if  obtained,  will 
apply  only  for  the  linearized  system.  This  is  considered  to  be  the 
best  the  software  package  can  do  in  the  circumstances.  The  problems 
inherent  in  simulating  an  unstable  system  on  the  analogue  portion, 
e.g.  limiting  because  of  particular  choices  of  scaling,  are  clearly 
problems  which  must  be  left  to  the  user. 


4.  DESCRIPTION  OF  THE  SOFTWARE  PACKAGE  AND  ITS  USAGE 

The  purpose  of  this  Section  is  to  provide,  with  reasonable  completeness,  the 
information  required  by  a potential  user  of  the  package. 

4.1  Description  of  the  software  package 

The  package  is  programmed  in  Fortran  for  an  IBM  370/168  digital  computer. 

It  uses  library  routines  from  the  Harwell  library,  for  finding  the  inverse  of 
a complex  matrix,  and  the  I.M.S.L.  library  for  computing  the  eigenvalues  and 
associated  eigenvectors  of  both  real  and  complex  matrices. 

The  program  is  best  described  with  the  aid  of  a flowchart,  shown  in 
figure  2.  The  matrices  Ai  , ...  D4  and  their  dimensions  are  first  defined 
and  printed.  The  eigenvalues  and  eigenvectors  of  Ai  are  computed,  and  the 
properties  of  the  real  part  of  the  eigenvalues  are  examined.  The 
eigenvalues  are  also  examined  to  establish  if  any  are  repeated.  Next,  the 
fundamental  matrix  for  Ai , and  the  matrices  G and  Q are  computed  for  a range 
of  step- lengths,  h,  chosen  by  the  user.  The  direct  or  numerical  integration, 
(N.I),  method  is  used  for  the  previous  step,  depending  on  user  choice,  or 
whether  the  eigenvalues  of  Ai  are  distinct,  (N.I.  method),  or  the  real  part 
of  any  eigenvalue  is  positive  (direct  method).  For  each  step-length,  h,  the 
partitioned  matrix  D is  computed,  and  its  eigenvalues  are  examined  to 
determine  the  stability  of  D.  The  main  program  prints  information  at 
intervals  throughout  each  stage  of  program  execution. 

The  Fortran  coding  and  a description  of  the  variables  used  are  presented 
in  Appendix  III  and  IV.  The  program  uses  double  precision  arithmetic 
The  matrices  A| , ...  A*,  D|  , ...  D* , G(h),  Q(h)  and  D(h)  are  real,  while 
the  eigenvalues  and  fundamental  matrix  for  ^ are,  in  general,  complex.  As 
a result,  some  of  the  subroutines  are  more  involved  than  might  otherwise  be 
expected. 
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4.1.1  Main  program 

The  main  program  accepts  the  input  data,  controls  the  main  branches 
of  the  analysis,  namely  which  method  is  used  to  calculate  G(h),  and 
Q(h),  and  prints  the  results.  It  makes  calls  directly  to  subroutines 
EIGEN,  FUND,  LAMBDA,  NUMERI,  DMATRX,  MODEIG,  and  OUTPUT,  which  are 
described  below. 

4.1.2  Subroutine  EIGEN 

This  routine  controls  the  computation  of  the  eigenvalues,  and 

the  matrix  of  eigenvectors,  C,  for  Ai  , via  subroutine  CNVRT. 

(C  = P*'  in  the  theory  of  reference  1,  Appendix  I).  EIGEN  calls  sub- 
routine DIST,  after  which  it  tests  if  the  numerical  integration 
technique  has  beer  chosen  by  the  user,  and  whether  X^  are  repeated. 

If  the  latter  situation  exists,  the  routine  forces  the  program  to  use 
the  N.I.  method.  If  neither  of  the  above  criteria  is  true  the 
routine  computes  the  inverse,  B,  of  the  matrix  of  eigenvectors,  C,  via 
subroutine  CNVRT,  in  preparation  for  computing  the  fundamental  matrix 
corresponding  to  Ai  by  the  eigenvalue  method.  Should  the  inverse 
routine  fail,  the  N.I.  method  is  chosen. 

4.1.3  Subroutine  DIST 

Specifies  the  real  pai i of  each  of  the  X.,  for  determining  the 
stability  of  G(h),  and  aiso  establishes  whether  the  X^  are  repeated, 
i . e . , whether  X . = X . for  i j . 

1 j 

4.1.4  Subroutine  FUND 

The  purpose  of  this  routine  is  to  define  the  fundamental  matrix, 

*(h),  by  direct  calculation,  using  eigenvalue  theory.  For  each 
value  of  h,  the  routine  computes  two  complex  diagonal  matrices 

diag  (exp  X^h,  ...  exp  Xj^h)  and  diag  (exp  — * “■  exp  i 

using  a complex  matrix  multiplication  routine,  DCMPY.  It  is  expected 
that  the  imaginary  elements  of  the  above  matrices  will  be  less  than 
lO"*  in  magnitude.  The  matrices  are  redefined  as  G(h)  and  Q(h).  A 
diagnostic  is  printed  if  the  imaginary  elements  of  the  complex  matrices 
do  not  conform  to  the  above  requirement,  but  the  program  continues  to 
execute . 

4.1.5  Subroutine  NUMERI 

This  routine  provides  the  control  for  the  numerical  integration 
routines.  The  numerical  integration  step-length  is  established  by 
routines  LAMBDA  called  from  MAIN,  and  STEP,  and  the  integration,  using 
the  "predictor-corrector"  method,  is  performed  in  subroutine  SINT. 

NUMERI  increments  time  up  to  the  required  value  of  h,  and  computes  G or 
Q depending  on  the  value  assigned  to  a variable  NAM. 

4.1.6  Subroutine  LAMBDA 

In  this  subroutine  the  analysis  of  Appendix  II  is  applied  to  establish 
the  nuoaerical  integration  step-length,  chosen  for  stability  and 
accuracy  of  the  solution.  Values  for  the  N.I.  step- length,  h,  are 
computed  as  a function  of  the  real  part  of  each  eigenvalue,  X^,  and 

must  also  meet  a criteria  related  to  the  imaginary  part  of  the 
eigenvalue  X^  (as  described  in  Appendix  II). 

If  this  criteria  is  met  for  all  the  eigenvalues,  the  smallest  h is 
chosen  for  the  N.I.  step- length. 
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If  the  above  criteria  is  not  met  for  one  or  more  of  the  eigenvalues, 
then  h|  is  chosen  for  the  respective  eigenvalues  according  to  the 
analysis  of  ^pendix  II,  and  the  N.I.  step-length  is  the  minimum  over 
all  the  eigenvalues  of  the  h and  h|'s  respectively.  h must  also  be 
at  least  as  small  as  the  minimum  interval  between  the  step-lengths 
chosen  by  the  software  package  user  for  examining  the  stability  of  the 
hybrid  solution. 

The  subroutine  is  entered  from  the  MAIN  program,  and  provides  the 


N.I.  step- length,  h 


STEP’ 


for  further  refining  in  subroutine  STEP. 


4.1.7  Subroutine  STEP 


Subroutine  LAMBDA  defines  the  N.I.  step-length  purely  on  the  basis 

of  an  accuracy  criteria,  and  the  largest  eigenvalue.  It  is  necessary 

to  establish  individual  N.I.  step-lengths  for  each  interval  of  stability 

step-lengths  chosen  by  the  user.  The  numerical  integration  step- 

length,  h , for  the  interval  (h. , h.  ,)  is  modified,  such  that 
STtr  1 1+1 


h 


S 


’'step 


± € 


where  e is  a small  quantity,  and  there  are  exactly  an  integral  number 
of  integration  steps  in  the  interval  (h^^,  Subroutine  STEP  is 

called  from  subroutine  NUMERI  whenever  a new  N.I.  step-length  is 
required.  The  subroutine  returns  the  integration  step-length  hg. 

4.1.8  Subroutine  SINT 

The  predictor-corrector  steps  for  numerical  integration  are  defined 
in  subroutine  SINT.  The  calling  vector  defines  IPASS  and  NAM,  to 
control  which  branch  of  the  routine  is’  to  be  executed,  i.e.  the 
"predict"  or  "correct"  step,  and  for  identifying  which  of  functions  G 
or  Q is  being  computed. 

4.1.9  Subrout ine  DMATRX 

The  partitioned  2H  x 2M  matrix  D is  defined  in  this  subroutine. 

It  calls  the  subroutine  D>WPY,  to  perform  multiplication  of  two  real 
matrices. 

4.1.10  Subroutine  MODEIG 

The  routine  controls  CNVRT  to  compute  the  eigenvalues  of  D,  and 
then  tests  that 


IPjl  < 1,  i » 1,  ...,  2M 

If  this  criteria  does  not  hold  for  all  p^,  the  routine  records  which 

of  the  eigenvalues  are  in  modulus  too  large,  and  their  respective 
values,  for  printing  by  the  main  program. 

4.1.11  Subroutine  CNVRT 

This  subroutine  controls  all  calls  to  the  library  routines.  It 
has  three  branches,  each  being  defined  by  a value  I,  assigned  in  the 
calling  subprogram  according  to  whether  eigenvalues  and  eigenvectors 
of  a matrix  are  required  (I  > 1),  the  inverse  of  a matrix  is  required 
(I  “ 2),  or  eigenvalues  only  are  required  (I  ■ 3).  EIGEN  calls  CNVRT 
for  I - 1,  2,  and  MODEIG  calls  CNVRT  for  I - 3. 


I 
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!'or  T ■ 1,  CNVRT  calls  RIGRP  from  the  IMSF.  library,  which  computes 
complex  eigenvalues  and  complex  eigenvectors  of  Ai  , and  a performance 
index  for  the  mathematical  accuracy  of  the  technique.  EIGRP  returns 
the  performance  index  and  CNVRT  prints  a diagnostic.  It  is  possible 
that  if  the  eigenvalues  are  not  distinct,  then  the  matrix  of  eigen- 
vectors is  incorrect.  The  software  package  prevents  this  matrix 
being  used  further,  by  forcing  the  program  to  use  the  numerical 
integration  method  for  computing  the  fundamental  matrix  4’(h), 
whenever  the  eigenvalues  are  repeated.  The  eigenvalues  are  not 
normalised  by  EIGRP. 

Por  I “ 2,  QWRT  computes  the  inverse  matrix  B of  the  complex 
matrix  of  eigenvectors  C,  using  the  function  MA23BD  of  the  Harwell 
library  routine  MA23AD.  The  library  routine  overwrites  the  input 
matrix,  and  returns  the  inverse  matrix  in  the  storage  provided  by  it. 

It  is  therefore  necessary  for  the  program  to  remember  C before 
calling  MA23BD. 

The  routine  MA23BD  makes  no  allowance  for  the  size  of  the  matrix 
in  the  dimension  or  common  statement  of  the  calling  subprogram  to 
be  bigger  than  the  matrix  row  dimension  in  the  calling  statement. 
Therefore  it  is  necessary  to  define  another  complex  matrix,  Y, 
having  a variable  dimension  in  the  calling  subprogram  CNVRT.  This 
branch  of  CNVRT  is  only  executed  with  IMETH  = 1. 

For  1=3  CNVRT  calls  EIGRP,  to  compute  the  complex  eigenvalues  of 
the  matrix  D.  This  call  to  CNVRT  is  made  from  subroutine  MODEIG 

4.1.12  Subroutine  D^•^PY  and  DCMPY 

Subroutine  DMIPY  performs  matrix  multiplication  of  two  real 
matrices.  The  matrices  need  not  be  square,  but  the  program  tests 
that  the  coltnn  dimension  of  the  first  matrix  equals  the  row  dimension 
of  the  second  matrix.  TTie  output  matrix  is  10  x 10,  even  though  the 
submatrix  containing  the  solution  is  x LLL  which  may  be  smaller 
than  10  X 10.  DCMPY  is  the  complex  version  of  DMMPY. 

4.1.13  Subroutine  (XITPUT 

This  subroutine  controls  the  printing  of  the  input  matrices 
Ai  , Di  , ...  D*  . 

4.1.14  Library  routines 

EIGRP  and  MA23BD  are  well  doctnaented  in  the  Computing  Services  Group 
Library  of  N.R.E.  The  Harwell  and  IMSL  libraries  are  both  extensively 
used  and  tested  by  other  computing  facilities,  and  have  performed  well 
in  this  software  package. 

4.2  Software  package  running 

This  Section  covers  the  data  requirements,  the  program  output,  and  the 

method  of  creating  a data  file  under  the  IBM  370/168  stq>plied  TSO  system. 

Operational  characteristics  are  explained,  and  the  identification  of  the 

W.R.E.  Library  version  of  the  software  package  is  supplied. 

4.2.1  Data  input  and  program  output 

Table  1 lists  the  input  parameters  and  the  form  in  irfiich  they  must 
be  presented,  but  further  discussion  of  the  parameters  is  presented 
b^low. 

Matrices  AI,  Dl,  ...  A4  must  be  defined  in  terms  of  their  sizes 
and  elements.  In  most  system  problems,  some  at  least,  of  these 
matrices  will  be  null,  and  the  others  will  very  likely  be  sparse. 

The  user  need  only  define  the  non-zero  elements,  and  the  respective 
row  and  coliani  of  each  element,  for  the  non-null  matrices. 
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TABLE  1.  DATA  REQUIREMENTS  FOR  SOFTWARE  PACKAGE 


Data 

record 


Paraaeter  value 


Parameter 

description 


Constraints  on 
parameters 


1 K,M  K > Row  size  of  A1  suitrix  K < 10  (nominally) 

M “ total  number  of  inte-  K ♦ 1 < M < 20  (nominally) 
grals  computed  analogue  If  no  integrals  are  assigned 
or  digitally  to  digital  computer  i.e. 

M = K,  then  M set  to  K ♦ 1 

2 NEL  Numbers  of  non -zero 

elements  in  A1 


3 


I,J,A1(I,J) , . . . I « row  number 
for  NEL  J = column  number 

elements  A1(I,J)  » non  zero 

element  of  A1  in  I^^ 

row  and  column. 


4 NEL 


Number  of  non- zero 
elements  for  D1 


I,J  are  integers 
A1(I,J)  is  a floating 
point  value. 

These  values  are  separated 
by  commas . 


. This  set  continues  until  Al,  Dl,  ...,  A4,  D4  have  all  been  defined.  If 
NEL  is  zero  for  any  of  these  matrices,  the  next  data  record  contains  NEL  for  the 
next  matrix  in  the  set.  The  order  of  input  of  the  matrices  must  run  Al,  Dl,  A2, 

D2,  A3,  D3,  A4,  D4.  The  data,  if  all  these  matrices  have  non-zero  elements,  are 
defined  in  records  2 to  17.  For  most  problems  there  will  be  less  records,  but 
NEL  must  be  defined  for  each  matrix,  even  if  it  is  assigned  the  value  zero. 


18  NT, EL, HINT 
(or  less) 


19  IMETH 
(or  less) 


NT  » number  of  user 
defined  step-lengths 
EL  • value  of  1st 
required  step-length 
HINT  - time  interval 
between  successive 
step- lengths. 


The  values  for  EL  and  HINT, 
are  required  if  the  step- 
lengths  are  at  equal  intervals. 
Otherwise  EL  ■ 0,  HINT  - 0. 

If  only  one  step-length  is 
to  be  defined  HINT  - 0. 


User  option  for  choos- 
ing direct  method  for 
computing  ^(h)  i.e. 
IMETH  • 1,  or  N.I. 
method  i.e.  IMETH  > 2. 


IMETH  ■ 1 may  be  over- 
written in  program  for 
N.I.  method. 

IMETH  > 2 is  considered  to 
be  accurate  enough,  in 
general,  and  is  recommended 
to  user. 


20  H(I),I-1..,NT 

(or  less) 


H(I)  are  the  user 
defined  step- lengths. 
NT  was  defined  earlier. 


This  record  is  (q>tional. 
It  will  not  be  needed  if 
H(I)  is  fully  defined  by 
data  record  18. 
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The  matrices  are  ordered  implicitly  Al,  Dl,  A2,  D2,  A3,  D3,  A4,  D4  in 
the  data. 

A variable  NEL  defines  the  number  of  non-zero  elements  for  each 
matrix  in  turn.  If  NEL  is  zero,  then  the  matrix  is  assumed  null,  and 
NEL  for  the  next  matrix  of  the  set  is  the  next  data  required.  When 
NEL  is  non-zero,  the  next  data  input  defines  all  the  non-zero  elements 
(NEL  of  them)  for  the  matrix  in  question. 

Consider  the  example,  in  which  there  are  4 non-zero  elements 


Al 


Assuming  the  siatrix  dimension,  i.e.  3 x 3,  has  already  been  defined, 
then  two  data  records  are  required  to  fully  define  Al.  The  first 
data  record  contains  the  value  of  NEL.  The  second  data  record 
contains  the  row,  column  and  value  of  each  non-zero  element  in  that 
strict  order,  although  the  order  in  which  each  element  group  is 
defined  is  unimportant. 

The  row  and  column  are  defined  by  integers,  and  the  element  value 
is  defined  by  a floating  point  constant.  The  two  data  records 
defining  Al  are: 

4 

1,1, 1.0, 2, 2, 2. 0,3, 2, 1.0, 3, 3,1. 


Following  the  matrix  definition  for  Al,  ...,  D4,  the  user  must 
define  a value,  or  values,  of  step-length  h,  for  which  he  wishes  to 
test  his  partitioning  scheme.  The  program  accepts  variables  to 
define  a range  of  values,  h,  at  equal  intervals,  or  a series  of 
values  h,  or  will  read  a single  value  for  h.  In  the  last  instance, 
the  software  package  provides  the  stability  for  the  single  value,  h, 
and  also  for  1.5  times  that  value.  This  may  give  the  user  a better 
idea  of  where  his  choice  lies  relative  to  stability.  If  the  step- 
lengths  are  at  equal  intervals,  the  package  computes  the  array  H(I), 
from  the  value  of  H(l),  and  the  time  interval  between  step-lengths, 
up  to  NT  values. 

The  user  has  a choice  as  to  which  method  is  used  to  compute  G(h)  and 
Q(h)  via  the  fundamental  matrix  4(h).  The  direct  calculation,  using 
eigenvalue  theory,  is  allocated  by  the  variable  IMETH  « 1.  The 
approximate  solution,  using  mmMrical  integration,  is  assigned  in  the 
package  by  IMETH  • 2.  Results  from  either  method  have  not  differed 
significantly  for  the  purpose  of  stability  analysis,  and  the  user  is 
advised  to  use  the  numerical  integration  method,  if  he  has  no  preference. 

Initially,  the  program  print-out  defines  the  input  matrices  Al,  Dl, 
...,  04,  together  with  the  values  of  h for  program  execution.  The 
output  then  defines  a performance  index  for  the  eigenvalue  routine, 
followed  by  the  eigenvalues  of  Al,  and  an  error  parameter  indicating 
success  (lER  ■ 0)  or  failure  (lER  > 0).  The  distinctness,  or  other- 
wise, of  the  eigenvalues  is  defined.  A comment  is  printed  concerning 
the  stability  of  the  matrix  G,  which  provides  the  user  a very  important 
guide  to  the  likely  success  or  failure  of  his  partitioning  scheme. 

This  stability  indication  is  independent  of  values  of  h.  Next,  the 
method,  by  which  the  matrices  G and  Q will  be  computed,  is  printed. 

The  preceding  inforMtion  is  printed  once  for  each  run. 
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The  next  block  of  output  is  repeated  for  each  value  of  h ascribed 
by  the  user.  The  step- length,  h,  in  seconds,  the  N.I.  step- length 
and  accuracy,  if  apiplicable,  and  the  matrices,  G,  Q and  D are  printed. 
The  value  of  the  error  parameter  for  computing  the  eigenvalues  of  D is 
given  and  the  eigenvalues  are  printed.  Finally,  there  is  a statement 
concerning  the  stability  of  the  scheme  for  the  present  value  of  h. 

4.2.2  Creating  a data  file 

In  order  to  run  the  software  package  under  TSO,  or  via  a batch  job, 
the  user  must  supply  a data  file,  and  link  this  file  with  the  Fortran 
program. 

For  most  user  problems,  the  matrices  Al,  D1  may  be  of  the  order 
6 X 6,  or  bigger.  The  record  length  required  for  specifying  the  non- 
zero elements  of  such  matrices  will  very  likely  exceed  80  characters. 
This  is  the  default  value  assigned  in  the  EDIT  coimand  or  by  the  JCL 
accompanying  the  data  cards , to  the  record  length  parameter  LRECL  (or 
LINE).  (The  8 characters  used  to  specify  the  line  number  are  included 
in  the  80) . A record  length  of  120  characters  is  considered  adequate 
for  problems  likely  to  be  encountered,  but  in  any  case  the  value  must 
not  exceed  255. 

The  BLKSIZE  (or  BLOCK)  parameter  of  the  EDIT  command  specifies  the 
maximum  length  (in  bytes)  for  all  the  records  of  the  data  set. 

Hence  the  BLKSIZE  parameter  must  be  an  integer  multiple  of  the  LRECL 
parameter.  The  default  parameter  is  3120,  (also  the  maximum  allowed), 
which  conveniently  is  a multiple  of  120,  but,  in  any  case,  a smaller 
number  may  be  specified.  For  the  test  data  sets,  described  in  this 
report,  the  BLKSIZE  parameter  has  been  specified  as  2640,  allowing  for 
a maximum  of  22  records  of  120  characters.  The  BLKSIZE  and  LRECL 
parameters  are  only  specified  when  a new  data  set  is  created. 

TIte  EDIT  command  for  creating  the  tracking  radar  simulation  data  set, 
SEEK4,  of  Section  5 is  given  as  an  example. 

EDIT  SEEK4  DATA  NEW  BLOCK(2640)  LINE (120) 


or 


EDIT  SEEK4  DATA  NEW  BLKSIZE (2640)  LRECL (120) 


4.2.3  Operational  characteristics 

This  Section  covers  the  general  operating  characteristics  of 
interest  to  a user. 

The  W.R.E.  Library  program  identification  for  the  software  package 
is  4154,  and  a user  will  be  supplied  with  the  data-set  name  of  the 
package,  and  any  other  JCL  routines  necessary  for  naming  the  program, 
on  request  from  Computing  Services  Groiq>. 

(a)  Running  time 

The  program  typically  coapiles  and  executes,  for  up  to 
twelve  values  of  step- length  h,  in  under  ten  seconds.  Exec- 
ution time  for  the  same  run  is  less  than  three  seconds. 

(b)  Diagnostics 

The  program  terminates  in  the  event  that  neither  the  N.I., 
nor  the  direct  method  for  coiputing  4>,  can  be  applied.  This 
will  occur  if  some  of  the  eigenvalues  are  repeated,  and  one  or 
more  have  positive  real  part. 

Most  of  the  printing  is  controlled  by  the  MAIN  program. 
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However,  there  arc  several  Instances  where  other  subroutines 
■ay  print  diagnostics.  In  sub-routine  FUND,  if  the  coaplex 
versions  of  the  matrices  G and  Q have  imaginary  elements 
larger  in  magnitude  than  10~’ , a diagnostic  is  printed.  The 
program  continues  executing. 

In  subroutines  DM4PY  and  DCMPY,  if  the  matrices  supplied 
are  not  coag>atible  for  multiplication,  then  a diagnostic  is 
printed,  and  the  program  terminates. 

Whenever  subroutine  CNVRT  uses  the  IMSL  routine  EIGRF  to 
compute  eigenvalues  and  eigenvectors  of  Ai  , a performance  index 
is  computed  defining  how  the  IHSL  routines  performed.  CNVRT 
prints  the  appropriate  diagnostic,  and  the  package  continues  to 
execute.  EIGRF  also  provides  an  error  parameter,  lER,  and  its 
value  is  printed  in  CNVRT.  lER  may  take  the  following  range 
of  values, 

lER  s 0 Routines  worked  normally. 

lER  > 128  The  numerical  solution  failed  to  converge 
for  one  of  the  eigenvalues. 

lER  < 128  One  of  the  parameters  in  the  calling 
statement  for  EIGRF  is  inconsistent. 

The  user  is  advised  to  regard  the  program  results  as  incorrect 
whenever  the  error  parameter  is  non-zero.  It  may  then  be 
necessary  for  the  user  to  refer  to  the  write-up  of  the  IMSL 
routine  EIGRF  in  order  to  satisfactorily  continue  with  the 
problem  solution. 

(c)  Replacing  software  package  routines 

Should  a user  need  to  modify  or  replace  any  of  the  routines 
of  the  software  package,  the  point  is  made  that  modifications 
or  new  routines  must  use  double  precision  arithawtic. 


5.  APPLICATIONS 


Any  physical  system  if  analysed  in  great  detail  is  non-linear.  Most  physical 
systems  are  also  time-varying,  but  if  changes  in  the  system  characteristics  are 
very  slow,  compared  with  variations  in  the  input,  a linear  time-varying  system 
can  be  approximated  by  a linear  time-invariant  one. 

It  has  been  assumed  here  that  a potential  software  package  user  has  derived  a 
linearized  version  for  the  system,  prior  to  considering  a hybrid  simulation 
model,  and  that  these  equations  are  in  the  form 


as  required  by  the  software  package. 

However,  it  sometimes  occurs,  that  in  neglecting  small  time  constants  in  the 
analogue  portion  of  the  model  certain  of  the  system  equations  are  algebraic 
equations  rather  than  differential  equations.  This  is  not  covered  by  the  theory 
presented  here.  Rather  than  extend  the  theory  to  cover  the  general  case  where 
equations  may  be  algebraic  or  differential,  it  is  considered  easier  to  reintroduce 
small  time  constants  for  the  purpose  of  stability  analysis. 

The  following  sections  deal  with  two  linear  systems  (see  reference  3 fbr  a 
detailed  derivation  of  the  linearization  process).  In  fact,  the  existing 
digital  simulation  model  for  example  2,  represents  the  control  fin  actuator,  for 
a missile,  as  responding__instantaneously,  to  produce  a fin  deflection,  i),  to  the 
autopilot  error  signal,  tf.  A lag  of  0.001  s is  introduced  to  express  the  system 
in  a form  suitable  for  software  package  implementation. 
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5.1  Tracking  radar  siaulation 

The  first  exaaple  describes  a (single  plane)  similation  of  a radar  guided 
weapon  hoaing  on  a target  (figure  3),  in  which  the  notion  of  the  nissile 
(under  autopilot  control)  is  simulated  on  the  analogue  computer,  while  the 
signal  processing  functions  of  the  missile  radar  receiver  and  the  dynamics  of 
the  radar  servo  and  pedestal  systems  are  represented  on  the  digital  coBq>uter. 
The  configuration  is  as  shown  in  figure  4. 

The  linearized  equations  which  describe  the  system  are 


'^'dm  “ ■ V 

T - -3(\(»p^  ♦ (3) 

i'j  * 7 (4) 

and 

r = G(S)  7 (6) 


where 


^01^  - dish-body  angle 


- missile  azimuth 


F(S) 


required  dish  rate  to  maintain  track 
missile  angular  velocity 
(1  ♦ Tjjj  S) 


1 ♦ T,,-  S ♦ T„,  S 


’D2 


D3 


and 


K(1  ♦ T S)(l  ♦ T..  S) 

■ («  • T„  S)(1  . T„  s S>) 

The  partitioning  scheme  is  such  that  and  7f  and  hence  are  evaluated 
and  integrated  digitally,  whilst  and  r are  evaluated  and  integrated  on  the 
amalofua  computer. 

It  is  necessary  to  reduce  the  second  order  equation 


F(S)  (♦j  - 
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to  two  first-order  equations.  Expanding  (2),  we  have 

'^DM  ^ ^D2  ^D3  ^DM  " ' '^M  ^D1  ’ ^D1 

Define 


z,  . - *I  • *„ 


Then  equation  (7)  reduces  to 


:!d2^ 


Integrating  the  above  equation,  we  have 

i , llsi^  ._1_Z,  -Zoi  ^ 

'^DM  T,3''dM  Tj,3"'  1^,3 

which,  together  with  equation  (8)  define  equation  (2)  as  two  first  order 
d.e.'s.  We  see  that  which  is  calculated  in  the  interval 

tj^  < t < tjj^j  on  the  digital  computer  is  a function  of  and  (^j)||(  from 

the  digital,  and  from  the  analogue  computer  in  the  interval  t„  < t < t„. 

M N"  1 N 

It  is  also  necessary  to  reduce  the  third  order  equation 

r = G(S)  7 


to  three  first  order  d.e.'s. 
Express  equation  (6)  as 


K(1  + S)  7 

1 ^ ^3  S 


’^A4  S * ^AS  S* 


and  define 


We  can  then  define 


K(1  ♦ T^j  S) 


and  hence 
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Equation  (6)  has  been  reduced  to  three  first  order  d.e.'s,  i,  Zi  and  Zj . 

Ne  shall  assiae  that  Zj  and  Zj  are  coaputed  on  the  analogue  coaputer,  and  can 

therefore  see  that  the  element  y is  computed  digitally  in  the  interval 

tjj  1 ^ passed  to  the  analogue  portion  via  the  D/A  interface  at 

time  tj^,  for  conputation  of  f,  Za  and  Zs  in  the  interval  tj^^  < t < 

Equations  (2)  to  (6)  can  now  be  expressed  by  a set  of  equations  defined 
in  (11). 


where 


Expressing  equations  (11)  as  d.e.'s  which  are  calculated  in  the  interval 

< t < t by  the  hybrid  coaputer,  we  have 
1 
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fH 

Q O 
H H 


{-•  >*>  “ 

< I M 
H X . < 

lO  K>  h- 


CM*  K) 

o a 
H 1h 


to 

< . < 
H H 


to 

< < 
H H 


<N  I/) 

h-'h** 


? < 

t-  H 


s r s 


< 

•o 


fH  K) 

0 a 
H H 

1 
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The  analogue  vector  x.»  and  the  digital  vector  are  defined  as 


^ ),  and  i = 
Zj 


It  remains  to  express  the  set  of  equations  defined  In  (12)  In  the  form 
defined  in  Appendix  I. 


0 ^A2  I ^A1 

■^AS  ^5 


- 1 - 


0 2i 


0 \Z3 


^A1  ^2 
^A3  ^A5 


0 0 0\ 


0 0 0 


0 0 0 


0 0 0 


^1  '^A2 
^A3  ^A5 
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i.e. 

ic  =>  Aijt  + 


and  Aj  a null  matrix.  The  equation  (13)  defines  the  analogue  portion  of 
the  partitioning  scheme. 

The  non-zero  elements  of  the  matrix  Dj  are  components  of  f,  Zi  and  Z3 , 
and  are  computed  on  the  digital  computer  as  shown  in  figure  4.  It  is 
immaterial  whether  the  elements  of  Dj  are  computed  digitally  or  on  the 
analogue  portion  (then  called  /^ ) since  the  two  matrices  occur  as  a sum 
(Aj  +02)  in  the  matrix,  D,  of  the  theory.  Similarly,  it  is  immaterial 

which  of  the  matrices  A3  or  D3  is  defined  in  the  expression  for 

Further,  expressing  (12)  in  the  form  of  Appendix  I for  ^ 


i 


0 0 0 

0 0 0 
0 0 Oj 


D3  ♦ D4  ^ (or  A3  + D4 


(14) 


with  A3,  A4  null  matriccK. 

Consider  as  an  example  under  consideration  the  situation  where 


D1  “ 

0.1 

D2 

0.22 

D3 

0.016 

K » 

3.52 
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= 0.023 


= 0.094 


Then  the  data  for  defining  the  configuration  for  the  software  package  is 

K = 4 

M = 7 

M-K  = 3 

A| , Di  are  4x4  matrices 

Aj , Dz  are  4x3  matrices 

Aj , Dj  are  3x4  matrices 


A* , are  3x3  matri  ces 


0 -5.8  11.4 

0 0 -43.5 


0 0 0 0 
-92.8  000 

459.1  000 

4.6  0 0 0. 


Aj  * fl 


-92.8 


-459.1 


Aj  a 5 


A*  ■ ® 


■13.8  6.3  -62.5 


Table  2 presents  the  data  in  the  form  accepted  by  the  package,  and 
Appendix  VI  shows  a sample  output  from  the  computer  run. 
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TABLE  2.  INPUT  DATA  FOR  EXAMPLE  1 

Parameter 
name 

K,M 

NEL 

A1 

NEL 
D1 

NEL  (A2) 

NEL 
D2 

NEL  (A3) 

NEL 
D3 

NEL  (A4) 

NEL 
D4 

NT,  EL,  HINT 
IMETH 

5.2  Height  control  system  simulation 

The  second  example  concerns  the  simulation  of  the  height  control  system  of 
a missile.  For  the  present  purposes,  the  system  is  assumed  to  have  the 
configuration  shown  in  figure  5.  In  this,  the  indicated  missile  height  and 
velocity  are  obtained  by  integration  of  the  missile  acceleration,  as  sensed 
by  the  inertial  platform.  Drifts  and  biases  are  corrected  by  comparison 
with  the  indicated  height  and  the  smoothed  output  of  the  missile's  altimeter. 
The  indicated  height,  Z,  is  then  compared  with  the  demanded  height,  Zj.,  and, 
with  the  vertical  velocity,  v,  is  used  to  form  the  demanded  acceleration, 

^ a^.  This  is  the  input  to  the  missile  autopilot  which  controls  the  missile 

aerodynamic  response  such  that  the  achieved  vertical  acceleration,  a^,  closely 
follows  the  demanded  input. 

One  approach  to  the  hybrid  solution  is  to  simulate  the  autopilot  and 
inertial  platform  and  a portion  of  the  aerodynamics  on  the  analogue  computer. 
Evaluation  of  the  missile  acceleration  and  the  derivative  of  pitch  rate  are 
computed  digitally,  since  they  are  functions  of  several  variables,  and 
integrated  on  the  analogue  portion  (figure  6).  The  above  configuration  is 
the  "text-book"  version  of  the  partitioning  problem  viz,  that  algebra  should 
be  solved  on  the  analogue  computer. 

The  analysis  of  reference  1 has  shown  the  potential  danger  of  applying  this 
"text-book"  approach,  if  the  resulting  scheme  evaluates  derivatives  on  one 
portion  of  the  computer  and  integrates  them  on  the  other.  This  "complete 
partitioning"  can  result  in  an  even  less  stable  division  than  evaluating  and 
integrating  derivatives  entirely  on  the  digital  computer,  user  Euler 
integration.  Nevertheless,  without  the  guidance  which  the  analysis  provides. 


Data  Record 

4,7 

7 

1.2. 1.0.  2. 2, -5. 8, 4, 2, 1.0, 2, 5, 11. 4, 3, 3, -43. 5, 
4, 3, -0.6, 2, 4, -10. 6 

3 

2,1,-92.8,3,1,-459.1,4,1,4.6 

0 

3 

2.1, -92.8,3,1,-459.1,4,1,4.6 
0 

3 

1.1,  -6. 3, 2, 1,-3. 0,3, 1,1.0 
0 

6 

1.1,  -13. 8, 2, 1,-3. 0,3, 1,1. 0,1, 2, 6. 3, 3, 2, 

-1.0, 1,3, -62. 5 

6.0. 038.0.001 
2 
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this  partitioning  scheme  might  well  appear  to  be  a sensible  one,  and,  as  such, 
provides  an  interesting  test  case  for  the  theory.  The  linearized  system 
equations  are  defined  below. 

Z * V 


I 

i 


V » 


a 


Z 


I 


a^  “ aio  >•  ai  Tj 

» aja  * T) 


d 


u 


♦ q 


17 


A4 

r 


(1 


‘ \2  «(>  • S)  ’ ■ “ 


i 


“d  “ “^AZ^^D  ' ^ 

For  stability  analysis  we  can  regard  Zp  as  zero.  Thus 

®D  “ ‘"^ZZ  ^ ■ "^AS  '' 


The  closed  loop  response  for  missile  height  and  acceleration  are, 

, 

^ “ (1  ♦ 0.85  S)(l  ♦ 0.15  S ♦ 0.039  S*) 


and 


3(1  - 0.0035  SI 

*Z  ■ (1  ♦ 0.18  S)(l  ♦ 0.016  S)  D 

In  the  above  equations  ai  , aj , aa  and  a*  are  functions  of  several  variables 
related  to  missile  incidence,  Hach  number,  and  air  pressure. 


‘d 


u 


‘^Al'-'^AS 

^Ar’'A2’^A3 


demanded  missile  acceleration 

demanded  missile  height 

missile  velocity 

gains  in  pitch  autopilot 

time  constants  in  body  rate  feedback  of 
autopi lot 
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Thus  the  linearized  equations  for  stability  analysis  are 


z » V 


V ■ aiO  ♦ 32 1} 


()  « ajo  ♦ a4»? 


a « - (ma  ♦ a2i})  ♦ q 


T I_(l  ♦ T S)(l  ♦ T S)  ^ 


K^sCa.a  ♦ a2t,)  - Ka2  " ' hi  ^ 


']-’r 


he  aust  reduce  the  third  order  d.e.  defined  by  equation  (19)  to  three  first 
order  d.e. 's. 

Put 


1 ♦ T..,  S 


P - - —0 

■ a • a 


A2  •A2 


he  will  assune  that  0 is  coaputed  on  the  analogue  computer, 
now  becoaes 


Equation  (19) 


-(1  > S) 


^4 

r |_C1  ♦ T 


- K (aio  ♦ aiq)  - K _ z - K. 


,J.i 


<■  ♦ '^>1 
(I  . T„  S) 


A3  *A3  ‘a3 


-8  ^A1  *A1  P r ’’^Al 

7 * T T * T — I ’ " T 

'a3  'A3  A2  'A3  L ‘a2  J 
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Substituting  for  8 in  equation  (20)>  now  have  a first  order  d.e.  for  f}. 


-r-  - •'as  ■ ''as  ■ 


''a2  * ■ "as 


Thus  all  the  equations  of  the  height  control  system  are  now  defined  as  a set 
of  first  order  differential  equations. 

The  partitioning  scheme  is  such  that  there  are  no  integrations  performed 
digitally.  At  time,  t^^  the  values  for  a and  t;  are  passed  to  the  digital 

computer  for  evaluation  of  (v)„  and  (4)vi  in  the  interval  t„  , < t < t„. 

These  derivatives  will  be  integrated  in  the  interval  t„  < t < t„  , on  the 

analogue  computer.  We  will  assume,  for  the  purposes  of  defining  A|,...,A4^ 


Di Da  that  the  coefficients,  ^ 

digitally  in  the  interval  tj^  ^ < t < tj^. 


ai  "aA  "as 


are  computed 


The  equations  defining  the  partitioning  scheme  are 


^ ^ ‘I 


a,  (a), 


T ^ " T 

‘a2  ‘a2 


"aa  "as 


1 _ ''a4''a5®* 


AS®*  "aA  ''a2  . "aA  "as  "aA  , 

— r ^ r-  ® 


V « ai  (a) 


♦ a2  (»»). 


^A1  ^A1  1 / ^A1 

T — T~  **  T~  V ' T — 
‘A3  ‘A2  'A3  \ 'A2 


■''A3® 


Ne  have 


q 

q and  ^ 
z 
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In  the  interval  t^^  < t < we  have 


fL 

u 

0 

0 

u 

0 

0 

i 

aj 

0 

0 

34 

0 

0 

0 

q 

0 

0 

0 

0 

0 

0 

0 

0 

■'^A4  '^AS 

0 

0 

■'^A4  "^AS 

0 

0 

0 

f| 

T 

V 

T 

*9 

(1 

0 

0 

0 

0 

0 

0 

z 

1 

0 

0 

32 

0 

0 

0 

K 

\ 0 

0 

0 

0 

0 

0 

0| 

i =■  Ai  + D,  j 


The  aatriccs  Aj , D2 , As,  D3 , A4  and  D4  are  all  null 
Consider  the  situation  where 
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Thus  the  data  for  defining  the  configuration  for  the  software  package  is 


Table  3 lists  the  input  data  for  the  software  package  and  Appendix  VII 
presents  coaputer  output  for  the  above  exaaple.  As  is  to  be  eiqpected  the 
Partitioning  scheae  is  shown  to  bo  far  froa  ideal,  and  the  results  show  that 
it  would  be  aatheaatical ly  unstable  for  a step- length  as  saall  as  12  as, 
even  though  the  bandwidth  of  the  acceleration  and  height  loops  are  0.9  Hz 
and  0.2  Hz  respectively.  The  general  experience  is  that  even  the  use  of 
prediction  will  not  greatly  help  in  these  sorts  of  application  which  ewloy 
"coaplete  partitioning". 


WRi;-TR-1«ft2(A) 


26 


TABLE  3.  INPUT  DATA  FOR  F.XAMI’LE  2 


Parameter 

value 

Data  Record 

K,H 

7.8 

NEL 

11 

A1 

1.2.1.0. 3.2.19.5.7.2.975.0.3.3,-0.5, 

7. 3. 75. 0. 4. 4, -1000. 0,4, 5, -151. 9, 4, 6, -151. 9, 

5.6.1.0. 4.7.70.0.7.7,-100.0 

NEL 

8 

D1 

1,1,-0.6,2,1,-58.4,4,1,1709.3,6,1,-193.8, 

1,4,-0.1,2,4,-61.6,4,4,304.3,6,4,-34.5 

NEL  (A2) 

0 

NEL  (D2) 

0 

NEL  (A3) 

0 

NEL  (D3) 

0 

NEL  (A4) 

0 

NEL  (D4) 

0 

NT,  EL,  HINT 

5,0.009,0.001 

IMETH 

2 

A better  partitioning  scheme  is  one  for  which  a and  are  evaluated  on 

z 

the  analogue  computert  and  the  inertial  platform  is  simulated  digitally. 

The  software  package  was  used  to  determine  the  stability  of  this  partitioning 
scheme  for  the  numerical  values  of  the  previous  example  (but  with  aj  » 0). 

It  was  predicted  by  the  software  package  that  instability  would  occur  for  a 
digital  computing  step-length  of  113  ms,  which  is  a considerable  improvement 
over  the  12  ms  of  the  first  scheme.  This  value  was  checked  against  the 
hybrid  simulation  by  incorporating  a variable  delay  in  the  digital  computer. 
For  the  same  values  it  was  determined  experimentally  that  instability  occurred 
for  a step- length  of  110  ms,  which  is  in  close  agreement  with  that  predicted 
by  the  software  package.  In  actual  fact  the  hybrid  simulation  was  inaccurate 
at  this  step-length,  even  though  stable,  highlighting  the  need  for  future 
work  to  consider  the  accuracy  of  hybrid  simulation. 

Finally,  by  using  a simplified  description  of  the  above  partitioning 
scheme,  the  question  of  stability  can  be  settled  by  direct  analysis,  as  in 
Section  5.1  of  reference  1.  For  values  of  K|  , Kj  6.6  and  T » 0.2,  and 
following  the  method  of  reference  1,  the  maximum  value  of  h for  which  the 
partitioning  scheme  is  stable  is  120  ms.  Thus  there  is  a good  measure  of 
agreement  between  all  approaches. 

6.  TEST  STATUS  OF  THE  SOFTWARE  PACKAGE 

The  package  was  extensively  tested  during  each  stage  of  its  production. 

However,  testing  the  complete  package  posed  a problem  since  a typical  system  for 
analysis  is  likely  to  be  quite  complex,  and  one  can  expect  the  D matrix  to  be  of 
the  order  14  x 14.  Finding  the  eigenvalues  of  this  matrix  by  hand  is  quite 
intractable. 

However,  it  is  possible,  by  reducing  the  time  constants  in  the  response 
functions  of  the  analogue  portion,  to  generate  'instantaneous*  responses  for 
that  part,  and  thus  produce  a system  which  can  be  analysed  theoretically.  One 
can  therefore  theoretically  analyse  a partitioning  scheme  for  such  a system, 
and  hence,  produce  a value  for  step-length  at  which  it  will  became  uistable. 
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Such  modified  system  equations  can  be  used  as  the  test  data  for  the  software 
package . 

The  homing  missile  simulation  of  the  first  example  was  chosen  as  a test  case. 
The  configuration  is  changed  slightly  with  the  dynamics  of  the  radar  servo  and 
pedestal  system  now  being  represented  on  the  analogue  computer.  The  configur- 
ation is  shown  in  figure  7. 

The  data  defined  in  Section  5.1  is  applied  with  the  time  constants  of  the 
transfer  functions  being  increased  by  an  order  of  10.  The  equations  for  the 
instantaneous  system  are 

Digital 


*i  ■ ■’('‘'dh  • V 

= 7 

Analogue 

(1  + 0.01 

'^DM 

1 ♦ 0.022  S ♦ 0.00016 

3.52(1  + 0.001  S)(l  + 

0.19  S)  7 

T ® 

(1  ♦ 0.0023  S)(l  + 0.054  S 

+ 0.00094  S^) 

M 

r 

These  equations  reduce  to  seven  first  order  d.e.'s  with  the  analogue  vector, 
being  uefined  by 


X = 


r 

Z| 

zj 

Z3 


and  the  digital  vector,  y,  being  defined  by  a single  element. 


(♦,) 


I 

I 

I 


! 


The  partitioned  system  can  now  be  defined  in 
for  the  software  package. 

In  vector  form  the  system  is  defined  for  the 


terms  of  the  matrices  suitable 
interval  tj^  < t < t^j^j 


J 
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I In  term  of  the  matrices  A .A4  .Di  , . . . ,D4  . 

f 

I £ - Ai  jc  ♦ Di  j ♦ A2 


with  Di  null,  and 

i ■ Dj  X N 


with  Aj.  A4,  Di  null. 

Solving  the  stability  problem  theoretically,  we  have  on  the  digital  computer 
in  the  Interval  t^j  < t < 


' i 
i 
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*r  - < * - <> 


using  Euler  integration. 

In  the  sa»e  interval  the  analogue  functions  are  con^iuted  so  fast  as  to  be 
considered  instantaneous.  We  can  therefore  represent  the  analogue  functions  by 
the  following  set  of  equations. 


r » -3.52  X 

» T 


At  time,  tj^^j 

> <)  (21) 

- 3.52  X (22) 

fron  equation  (23),  we  have 

< = C * <»> 

Thus 

- 3h  (25) 

and 

^ - 3h(^<^'^)  (26) 

Equations  (25)  and  (26)  are  recurrence  relations,  which  together  define  the 
coaplete  systea.  It  is  clear  that  the  stability  of  equation  (26)  depends  on 
the  stability  of  equation  (25),  i.e.  stability  is  determined  froa 
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Clearly  a solution  to  equation  (25)  is 

^ = a,  0^*  > 32  02  (27) 

where  0i  and  02  are  the  roots  of 

0*  - 0 + 3h  = 0 (28) 

i.e. 

„ 1 ± n/  1 - 12h 

01,2  = ^ 

If  h < Yj,  the  roots  are  real,  with 

101,21  < 1 

Hence  equation  (27)  is  a stable  solution  to  equation  (25) . 

If  h > the  roots  are  complex  conjugates,  and  for  stability  of  equation  (27), 
their  product  must  be  in  modulus  less  than  unity,  i.e.  from  equation  (28)  we  need 

3h  < 1 

which  means,  for  stability  we  require 


The  software  package  was  run,  and  predicted  the  system  would  go  unstable  at 
exactly 


h = 0.33  s 


7.  CONCLUSIONS 

The  analysis  of  reference  1 has  shown  that  the  problem  of  stability  of  a given 
hybrid  simulation  can  be  solved  in  advance  of  programming,  provided  one  can  find 
a linear  approximation  to  the  response  of  the  system  under  consideration.  For 
large*sized  problem.s,  however,  the  steps  im  Ived  in  applying  the  direct  analysis 
may  be  intractable. 

This  report  describes  a software  package  which  implements  the  stability 
analysis  procedure,  thus  providing  hybrid  cong>uter  users  a quick  and  flexible 
means  of  determining  in  advance  the  stability  of  a proposed  partitioning  scheme. 
In  addition,  the  report  attempts  to  provide  guidelines  for  interpreting  results 
obtained  by  the  software  package. 


- 31  - WRi:-TR- 1882(A) 


8.  ACKNOWLEDGEMENT 

The  author  would  like  to  thank  Mr  M.J.  Fitzpatrick  of  Electronic  Warfare 
i • Studies  Group  for  his  assistance  in  the  programming  task. 

i 


WPR-TR-1882(A) 


32 
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A|  , Aj  , Aj  , A4 
B 

C 

D 


D|  . D7  . Dj  , D4 

FfS)! 

G(S)j 

G 

H(S) 

'k 

'^’'^Ar*'A2  ] 

'^A3*'^A4*'^A5  * 

N 

o 

P(S) 

Q 

^Al'^A2’^A3’] 
^A4’’*^A5’  [ 

^D1’^D2’^D3  J 

X(t) 


Z|  . Zj  , Zi 


“Z 

®l  « ®1  » *3  t ®4 


h.  h,. 


h 

o 


i 


NOTATION 

MXM  matrix  describing  linear  system  response 
matrices  used  in  partitioning  of  A 

matrix  used  in  partitioning  of  A,  also  inverse  matrix  of 
eigenvectors  of  Ai 

matrix  used  in  partitioning  of  A|  also  matrix  of  eigen- 
vectors of  Ai 

stability  matrix 

matrices  used  in  partitioning  of  A 
transfer  functions  used  in  example 

response  matrix  for  system  programmed  on  analogue  portion 
transfer  function  used  in  example 
unit  matrix  of  order  K 

gains  used  in  examples 

number  of  integration  steps  in  interval  (o,h) 

transfer  function  used  in  exan^le 

response  of  the  analogue  system  to  unit  step  input 

time  constants  used  in  examples 

unique  fundamental  matrix 
missile  height 
demanded  missile  height 

used  in  example  to  reduce  n^^  order  d.e.  to  n first 
order  d.e. 's 

maximum  element  in  matrix  A^ 
demanded  missile  acceleration 
missile  acceleration 

functions  of  several  variables  related  to  missile  incidence, 
Mach  number  and  airpressure 

hybrid  simulation  step-length 


#X 
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*'S’^STEP 

q 

r 


N 


u 


net) 

V 


a 

0 


0X.01 

y 
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numerical  integration  step-length 
missile  pitch  rate 
missile  yaw  rate 
time 

missile  velocity 
system  input 

missile  vertical  velocity 

vectors  used  in  describing  system  state 

error  function  for  fundamental  matrix 

null  matrix  of  order  K 

fundamental  matrix 
At 

equivalent  to  e 

pitch  incidence  in  example 

used  in  example  to  reduce  n^*'  order  to  d.e.  to  n first 
order  d.e. 's 

zeros  of  stability  quadratic  in  example 
receiver  output  in  example 

used  in  example  to  reduce  n^^  order  d.e.  to  n first  order 
d.e. 's 

accuracy  of  numerical  integration  solution 
pitch  wing  angle  in  example 
pitch  autopilot  error  signal  in  example 
eigenvalue  of  partitioning  matrix  Ai 
eigenvalue  of  stability  matrix  0 

time  nerements  for  numerical  integration  interval  (o,h) 
radar  dish  angle  in  example 
integral  of  receiver  output  in  example 
missile  attitude 
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APPENDIX  I 

ANALYSIS  OF  THE  GENERAL  PARTITIONING  PROBLEM 


Consider  the  problem  of  solving,  by  hybrid  computation,  the  system  of 
differential  equations 


i * + Cii(t)  (I.l) 

where  is  a vector  of  order  M,  and  ii(t)  is  a vector  function  of  time.  This  set 
may  represent  the  whole,  or  only  part  of  the  system  being  simulated,  but  it  is 
sufficient  in  the  present  context  for  investigating  the  stability  of  integration 
processes . 

For  the  purposes  of  stability  analysis,  one  can  disregard  the  forcing  function, 
and  set 


ii(t)  = £ 

in  equation  (I.l),  and  henceforth  be  concerned  with  the  system 


i = B*  (1.2) 

The  general  solution  to  solving  equation  (1.2)  on  a hybrid  computer  is  to 
allow  for  integrations  and  evaluation  of  derivatives  to  be  performed  digitally 
or  on  the  analogue  computer  without  restriction.  The  range  of  possibilities, 
then,  is  as  shown  in  figure  1. 

If  K integrations  are  carried  out  on  the  analogue  computer  and  M-K  on  the 
digital  portion,  then  equation  (1.2)  can  be  re-written  in  the  form 

i.  = A,y.  ♦ ♦ D,y.  + Dji 

i = ^3)L  * A4i  + DaX.  + D4i  (1.3) 

where 


and  the  vectors  A|)r.  Aj^,  A3X,,  Kti,  are  computed  on  the  analogue  portion  while 
DiK,  D)X  ^nd  04^  are  determined  digitally.  The  matrix  B as  defined  in 

equation  (1.2)  is  then  assumed  partitioned  such  that 

B = A ♦ D, 

where  A,  D are  matrices  comprising  those  elements  of  B which  are  evaluated  on 
the  analogue  and  digital  computers  respectively.  Further  A and  D are  assumed 
partitioned  such  that 
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D|  D, 

D3  D4 


For  the  purpose  of  this  analysis,  the  digital  to  analogue,  and  analogue  to 
digital  interfaces  are  assumed  to  be  sampled  simultaneously. 

Assume  4 is  sampled  by  the  analogue  to  digital  converter  at  time,  t^^,  and  the 

digital  computer  step-length  for  each  cycle  is  h seconds.  Then  at  time 

the  quantity  is  fed  through  the  digital  to  analogue  interface  for  determining 

a in  the  interval  tj^^^  < t < 

Euler  integration  on  the  digital  computer  gives  from  (1.3) 


Vl  = -^N  ^ h(A3X^  ^ A4i^_j  ♦ Day^  * D^^^) 


(1.4) 


In  the  same  interval,  tj^  < t < i.  can  be  integrated  on  the  analogue  by 

solving 


i = A,)r  ♦ A2i^  + D.y^  j + 

In  C((uation  (f.i;)  Aiy,  varies  with  time,  while 

is  constant  over  the  interval  tj^  < t < tj^^j. 

The  solution  to  equation  (1.5)  is  given  by  (see  Appendix  V) 


(1.5) 


= 4.(h)  ‘h-' (o)  y^ 


h 

* f «>>({)  .1.-' 

•>  o 


(O)  di(A2ij^  + (1-6) 


which  reduces  to 


where 


G » «l>(h)  <!>■’  (o) 


G(n  dt 


(1.7) 


and  4>(t)  is  a fundamental  matrix  for  Ai  . 
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Equations  (1.4)  and  (1.7)  can  be  cond)ined  into 


where 


Setting 


QD, 

hA4/ 

U-i 

equation  (1.8)  becomes 


G Q(A2  + Ih) 

h(A3  ♦ D3)  + hD^ 


/qd,  e,  \ 

This  system  can  be  further  transformed  by  setting 

% ■ ("^1 


* % 


and  then  equation  (1.9)  reduces  to 

/ G Q(A2  + Dj) 

/ h(A3  D3)  I^_j^  + hD* 


\ ^ 

\ 


0 C|(A,  . 0.  )\  /)[^\  /qo,  /xn., 

h(A.  . D.)  . hD.  ) Lj  \ ^ hA.]  \i„., 


^ is  a null  matrix  of  order  Kx(M-K) 
^ is  a null  matrix  of  order  (M-K)xK 


(1.8) 


(1.9) 


(I. 10) 


(I. 11) 
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NUMERICAL  INTEGRATION  STEP- LENGTH  FOR  DETERMINING  4>(h) 


It  is  required  to  deteraiine  the  solution  to  the  fundamental  matrix  equation 
using  niflierical  integration  over  the  interval  (0,h  ),  i.e.. 


where 


The  accuracy  and  stability  of  the  nmerical  solution  depend  on  the  numerical 
integration  step-length,  h,  which  will  be  determined  from  the  eigenvalues  of  B. 
The  numerical  integration  method  chosen  is  the  predictor- corrector  method* 


For  the  predictor- corrector  method  at  the  (N+l)  step  we  have 
Predict  pass: 


If  the  eigenvalues  of  B are  distinct,  then  we  can  find  a matrix  P,  such  that 


where  are  the  eigenvalues  of  B. 

(If  the  are  not  all  distinct  we  need  only  consider  the  reduced  set  of 

recurrence  relations  for  which  the  eigenvalues  are  distinct) . 

Setting 


we  obtain 


For  the  k column  we  have  for  the  predict  pass 
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Correct  pass: 


Vi  ■ 'S.*I<\*Vi) 


The  question  of  stability  can  be  determined  by  analysis  of  the  stability  of  the 
set  of  recurrence  relations 


Putting 


we  have 


(II. 1) 


and  require 


K h' 

0 = I - *-\- 


101  < 1 


for  a stable  solution. 

Hence  the  stability  boundary  is  defined  by 


K 

f(h)l  = I 1 + Xj^h  I = 1,  k 


= 1,...,  M 


(II. 2) 


and  is  shown  in  figure  8. 

It  is  necessary  for  the  accuracy  and  stability  required  in  the  present 
application  for  h to  be  chosen  so  that  X^^h  lies  well  inside  the  stability  boundary 

for  all  k. 

Consider  the  real  part  of  the  eigenvalues,  i.e.. 


Then  we  have  stability  if. 


u = Re(Xj^) 


uh  < -2 


In  order  to  achieve  uh  well  inside  the  stability  boundary  we  choose  h so  that 
we  have 


uh  ■ -0.2 


(II. 5) 


41 


WRE-TR-1882(A) 


To  further  ensure  accuracy  and  stability  it  is  necessary  that  the  above  value  < 
of  h satisfies 

j 

! 

where  i 

I 

V = Im  (Xj^)  ^ 

! 

w_(h)  = ^\/-8uh  (see  below)  , 

•»  i 

If  equation  (1 1.4)  is  not  satisfied,  then  it  is  necessary  to  define  a new  value 
of  numerical  integration  step-length,  hi , such  that  hi  satisfies  the  accuracy  i 

and  stability  criteria.  We  require  an  h|  such  that  uh|  are  small,  and  vht  lie 
well  inside  the  stability  boundary  for  all  the  eigenvalues  of  B.  [ 

Putting  * 

j 

X,  = u + iv 
k 

in  equation  (1 1. 2),  we  have  for  the  stability  contour  j 

f(h)  = 1 + uh  ♦ + i vh(l  + uh)  j 


For  complex  eigenvalues  we  are  interested  in  those  for  which  uh  is  small. 
Thus  equation  (I I. 5)  reduces  to 


2uh  ♦ v*h* 


uh  ♦ 


0 


Setting 


a ■ idi 
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wo  have 


If 


2a  + v'h^a  ♦ -5^  = 0 

j -4a  ± y/l6a^  - 32a 
(vh)"  = — 2 

= 2l  a|  ± v/8|a|  + 4a^ 
vh  = v/2l  a|  ± v^l  a|  + 4a^ 


Hence  for  uh  small,  vh  is  also  small,  and  in  equation  (I I. 5)  we  can  ignore  the 

term  in  v*h^,  but  not  the  4^^  order  term  since  it  is  independent  of  uh.  Thus 
the  stability  boundary  in  equation  (II.S)  can  be  expressed  by  the  relationship 


2uh 


-jvlh! 

4 


Suppose  we  choose  a step- length  hi , such  that 


w = vhi 


On  the  boundary  we  have 


Wg  = V-8uhi 


We  require  the  step-length  to  be  chosen  so  that  w is  well  inside  the  boundary, 
thus  we  choose  h|  such  that 


i . e . , we  have 


Wg(h| ) 


vh| 


< 0.2 


= 0.2 


^\/-8uh| 

Hence,  h|  is  chosen  by  application  of  the  expression 

hi  - X (0.2)’ 
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f 


t 


I. 


LISTING  OF  SOFTWARE  PACKAGE 
IMPLICIT  ReAL<'‘<(A-H),KEAL«0IO-Z) 

CUMPLEX«16  L( in),C( lOflUi tD( 10«10)«FI ( 10tlO)«PHll iO*lO) «GGI 
lQg< 10*101 •0FI(10tlQlfF(20l«NW(l53) 

COMMON  E«C*u*FI tPHl »GG*UO«UFI *F 

COMMON  A 1(  10*lo)fA2(  10*10)* A3(  10*10)  *A<i(  10*10)  *0i(  10*10)  « 

1021 10*10) *03 ( 10*10) «04( 10*10 )*H1 20) *901 t 10*10 )*wD2 ( 10*10) 

2*All(10*l0) 

COMMON  H(200)*2Z(10*ia)*Ll(lU)*EK(10)*G(10*l0)*C(10*10) 

COMMON  0«  20,20)  * T 1 M t * riS  *HS T bP* FPS  *U 
COMMON  PFI  I 10*  I0)*uri  ( 10*10)  *XINTA1(  10*  10  ) *X  1NTA<:  I 10*10) 
i*YFl( 10*10) ,0YF1( 10*lO)*YlNT Alt  10*10) *YINTA2( 10*10) 

COMMON  K*M*MMa*  UIG* INJEXl 10) * JCOUNTI 10) 

1*KCUUNT*IMFTH*  nt*mtwo 

DIMENSION  I IK.(  10)  * 1K(  10)  *OUM(  100)*X(  10*10  ) 

C read  dimension  of  MATRICES  DEFINING  THE  PARTITIONING  SCHEME 
READ! S*«) K*M 
MMK=M-K 
MTWU=M*M 

C ZERO  matrices  Al*...*0^. 

00  3 I=1*K 
nu  2 J=1*K 
All 1*J)=0. 

2 01(l*J)s0. 

00  3 J=1*MMK 
A2(  1*0) =0. 

D2( I* J) :0. 

A3( J*  1) >0. 

3 n3<J*I)-0. 

nu  <.  I=1*MMK 
00  M J>1*MMK 
A<*ll*3)sO. 

A 0A(1*J)*0. 

C READ  NON-ZERO  ELEMENTS  OF  Al* *0A. 

REAO(S*«)  NEC 
IFINEL*E9.0)G0  TO  5 

READ15*«) (I1K(L)*1K(L)*a1 (I1K(L)*1K(LI ) *L=1*NEL) 

5 READIS*«)  NEL 

1FINEL.E9.0)G0  TO  6 

READ! S*«l(IIK(LI*IK(L)*Ol<IlK(LI*IK(L)) *1=1* NED 
b REAOIS*«)  NEC 

IFINEL.EO.OGO  TO  7 

RtAD(3*«) (IIKIL)*IKIL)*A2(I1K(L)*1K(L) )*L=1*NEL) 

7 REAUf3*«)  NEL 

IFINEL.EO.OIGQ  TO  8 | 

READ! 5*0) (I1K(L)*1K(L)*02(IIK(L)*IK(L)) *L=l*NtL)  ! 

8 REACI5*«)  NEL  I 

IF(NEL.E0*0)GU  TO  9 I 

READIS*»)(11K(L),IK(L)*A3II1K(L)*1K(L)) *L=I*NEL ) i 

9 READI3*»)  NEL  | 

IF(NEL.EQ.O)GO  TO  11  < 

REAQ(5*«) (IIK<L)*IK(L)*D3(IIK(L)*IK(L) ) *L=1*NEL) 

11  REAU<S*«)  NEL 

lFINEL*tQ*0)GU  TO  12 

READIS*«) 1I1K(L)*IK(L)*AAIIIK(L)*1K(L) ) *L=I*NEL) 

12  REA0I5«»)  NEL 

IFINEL.EQ.OIGU  TO  13 

REA0(5««) niK(L)*lKIL)*DA( I1K(L)*1K(L) ) *L*I*NEL) 

13  REA0I3**)  NT*ELfHlNT 
REAU(3*«) IMETH 
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C tOMPUTt  UR  RfcAU  THt  KANGfc  UF  H FUR  TlSTING  STAbILITY. 
DU  10  I=ltNT 

10  HI1)=EL  ♦ DFLUATI  I-i)<=HlNT 

IFCHCND.EQ.OjRfcAutStO)  ( H ( I ) ♦ i = 1 , NT  ) 

C StT  UP  STEP  lengths 

IFINT.NE.DGU  TU  18 
2 » =1.5D0*H(  1) 

NT  = 2 


C 


WRITE  UUT  INPUT  DATA 

lb  WKlTE(bfl9) 

19  FURMATl 2Xt* INPUT  UATA*,/,tX 
^RI TE (6*1 S)KtM 
WRITE((>t21)KtK 

21  FUKMAT(2X,»a1  MATRIX-I • t I2» 
CALL  OUTPUT  (Al*K»KfX) 
WRIT£(6t22»K,K 

22  FURMATl2Xf *01  M ATKI X- I • . I 2 » 
CALL  OUTPUT  (Jl,KfX»X) 

TE  (6t23)KtMMK 

2J  FuRMAT(2X,'A2  MATRI X- I • 1 1 2 t 
CALL  OUTPUT  lA2»K,HMKfXI 
WKnE(6,29)Kf 

29  FURMAK  2X,'u2  M ATRI  X- I • » I 2 t 
Call  output  (D2»k,mmk,x) 
WKlTEl5,23)MMKf K 

2S  FORMAT! 2X.'A3  M ATR I X- I • ♦ I 2 t 
CALL  OUTPUT  (A3»MMKfK,X» 

WRI  TEIf>«26»MMKf  X 

20  FuRMATI 2X »»U3  MaTKI X- ( • » I 2» 
CALL  OUTPUT  IU3«MMK,KtX) 

WRI TO (6f2  TTMMKfMMK 

27  FORMATl 2Xf • AA  MATRI X- ( • ♦ I 2f 
CALL  OUTPUT  (A4*MMK,MMKtX) 
WRI TE (6«2d)MMK,MMK 

2d  FURMArC2X,*04  MATRIX-l',I2t 
Call  output  (uAtMMKfMMKfX) 
WRTTE(6t20)NTt(H(ntI  = ItNTI 


10« IH. J ,/) 

X* ♦ I2t • ) • ) 

X» , I2» • I • ) 

X*  »I2t'  » • ) 

X*  ,l2tM'  ) 

X*  tl2»'  ) • I 

X*  tl2t‘  I • » 

X*  »I2t  • I * I 

X*  f I2» • ) • ) 

ARt-'*/»in(2XtF6.A) ) 


lt>  F0RMAT!2Xt  •K  = ' t I 2 t 2X  , • M=  • » l 2 ) 

20  FORMAT  12X,*TH(  •.I2»*  STEPlEnGTHS 

WRI TE(6tl02) 

WRITEI6,l<i)  IMtTH 

19  forhat(2x,»ih£Th=»»  in 


UU  16  I=1«K 
00  16  J=nK 
16  All  I I»J)  ::A1  ( n J ) 

c compute  the  eigenvalues  anu  matrix  of  eigenvectors  fop  A1 

WRI  TC (6  ♦ 1 71 
17  FORMAT! IH  I » 

CALL  EIgEN 
WRI IEI6f 102) 

102  FuRMAT!1II  ) 

wR1Te(6« lOS) 

105  FOkMATIIOX.'EIgFnVALUES  of  A1  ICQMPLEX  EL EME NT S) • / t 1 OX * 35(lri-)) 
WRlTE!6t llul ir ! I )« I=1«K) 
no  format  (9  !2XfM  •»Lll.9f*»’»cll.9»')*)f/) 

WRITE!6«125) 

WRirEI6»lU2) 

C TEST  THE  ZEROS  OF  AltlE  WHEThEk  THE  E/VALUeS  HaVE 
C real  PART  NON-POSITIVE 
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rsT=o 
IMST=0 
INST=0 
00  36  I=1,K 

IF(  JCUU.NT  ( I l-l>  Jlf  33*  35 
31  lST=IbT*l 
GO  TO  36 
33  IMST=IMST*l 
GO  TO  36 

35  INST=INST»1 

36  CONTINUE 

IK  INST.NE. JJGO  Tu  60 
IF( IMiT.Nt.OIGO  To  50 
WRITEl6tA0) 

‘tO  F0RMAT12X*  6<*HZtRUS  Ot-  AT  HAVE  -VE  REAL  PART.O  MATRIX  WILL  bE  STAb 
ILE  FOR  ALL  H ) 

WRlTE(of  102) 

GO  TO  70 

50  WkITE(6»55)  IMST 

55  FORMATJ2X,I2*74H  OF  THE  ZEROS  OF  A1  HAVE  ZERO  REAL  PART.G  MATRIX  W 
liLL  BE  MARGINALLY  STABLE  ) 

WRITE! 6* 102) 

GO  TO  70 

60  WKITE!6*65)  INST.lMiT 

65  F0KMAT(2X,  VHTMERE  ARE*I3*29H  ZEROS  WITH  ♦VE  REAL  P AR T * AND* 1 3 » 56HZ 
lEROS  WITH  ZERO  REAL  PAPT.THOS  G MATRIX  WILL  RE  UNSTABLE  ) 
nR1TE(6* 102) 

IMETH=l 

IF( O.GE.O.UO)GO  TU  70 
WRITF(6,82) 

STOP 

70  CONTINUE 

C FOR  DISTINCT  EIGENVALUES  THE  FUNDAMENTAL  MATRIX 
C IS  COMPUTED  USING  THE  MATkIX  OF  EIGENVECTORS  OF  AI* 

C OR  USING  NUMERICAL  INTEGRATION-  USFK  CHOICE. 

C FOR  REPEATED  EIGENVALUES  THE  FUNDAMENTAL  MATRIX 
C IS  COMPUTED  USING  NUMcRICaL  INTEGRATION  ONLY 
IF!  lElG.EO.OCU  TO  7l 
wR1TE(6*75) 

75  F0RMATI2X*35HTHE  EIGENVALUES  OF  Al  ARE  REPEATED 

I*G  AND  U WILL  RE  OBTAINED  BY  NUMERICAL  INTEGRATION*  ) 

WRITEIb* 102) 

IF( INsT.LE.U)gO  Tu  72 
WRITE(6*82) 

82  FORMAU2X**aOTH  N.I  AND  DIRECT  METHOD  FUR  COMPUTING  FUNDAMENTAL* 
I**  MATRIX  WILL  FAIL*) 

STOP 

71  WRITE(6*73) 

73  FURMAT(2X** THE  EIGENVALUES  OF  Al  ARE  DISTINCT*) 

WRI rE(6«102) 

72  LL*NT 

1F( IMeTH.EO. 1 )CU  TO  74 
CALL  LAMBUA 

74  CONTINUE 

IPASS>0 

DO  100  INsl*LL 
TlME«Ht IN ) 

dO  IF( IMETH.EQ.2)CU  TO  84 

C COMPUTE  THE  FUNDAMENTAL  MATRIX  FOR  A1*AND  THE  MATRICES 
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C o AND  u F(JK  THl  KANGfc  UF  H SPECIFIED  AS  INPUT 
CALL  FUND 
WKI TE (6f I02> 

MKlTE(6tIll<»)TINC 
lli<»  KJKMAT  ( 3X» 'Hs*  ,Fc>.<.f  / » 

GO  TO  8t> 

C COMPUTE  GIH),0(H»  OY  NUMERICAL  INTEGRATION 
6^  CALL  NUMERI ( iPASStINtHiNT) 

WKI TE(6tl02) 
hRITEI6,93»TIME,HS 

93  FURMAI(2Xt*Hs*,F6.Af • NUMfcKICAL  INTEGRATION  STEPLENGTH  IS'fF9,6/) 
86  WKITE(6»11S) 

115  FOKMATUOXf  ‘oMATRlX*) 

OU  132  II=liK 

WRITE! 6t ill ) II t (b( II t Ji t J= i«K) 

III  FORMAT! 1H0»2X» 'ROW  * t I 2 t 10 ! I X, E 1 1 . A H 

132  continue 

wRITE!6t 1251 
WRITE! bf  112) 

H2  formaT!iox»  'u  matrix*) 

DO  133  n = l«K 

WRITE! 6t 11 1 ) II t ! U!  1 1 « J) t J= I tK) 

133  CONTINUE 

125  FOKMaTIIH  tl20llH-)) 
wRlTE!b«125) 

C COMPUTE  THE  MATRIX  0 OF  ORDER  2M 
CALL  DMATRX 
wRITE!6t 130) 

130  FORMAT ! lOX,8Hn  MaTkIa  » 

DO  131  lI=ltMTWU 

WRITE!6tll3) I1«!D! II*J)«JslvMTWU) 

113  FORMAT!  lH0f2X»*KUW  ' 1 1 2 « 1 0 ! I X t E 1 1 . <r ) «/ t 9X,  I U ! IX  * E 1 1 . A ) ) 

131  CONTINUE 

mRITE!&» 12?) 

135  F0KMAT!IM  fl20!lH«)) 

C TEST  THAT  THE  EIGENVALUES  UF  0 ARE  IN  MODULUS  LESS  THAN  UNITY 
CALL  MODEIG 
nRlTE!6«lAU) 

140  FOkMaTHOX, 'EIGENVALUES  OF  D !COMPLEX  ELEMENTS  )•  t/ ♦ 

110X,35!  IH-) ) 

WRITE!6t 110) !F I II) «I I = 1*MTwO) 

IF!KCOUNT.GI.O)GO  TO  67 
WKlTE16fR5)TIME 

65  F0RMATI2X,  •ALc  THE  EIGENVALUES  UF  0 HAVE  MODULUS  LESS  THAN  OR  Fj 
lUAL  TO  UNITY  FOR  H=**F6.4f  • SCHEME  IS  STaBlE  FOR  CHOICE  OF  H*) 
WRITE!6»I35) 

GO  TO  lUO 

87  WRITE !6*90)TIHE tKCOUNT 

90  FuRMAT!2Xt  bHFOR  H= F6 .4 » iHt 1 3 ♦ • OF  THE  EIGENVALUES  UF  D HAVE  MOD 

lUcUS  GREATER  THAN  UNI TY. THEREFORE  CHOICE  IS  UNSTABLE  FOR  THIS  H* ) 
WRITElbtl 10) !F! INuEXi II ) ) 1 1 1 ^ 1 1 KCDUNT ) 

WRI TE !o*l 35 ) 

IDO  CONTINUE 
STOP 
END 

SUBROUTINE  EIGEN 

IMPLICIT  REALP8! A-M) ♦ realms !0-Z ) 

COMPLEX*! 6 E( 101 tC( 10,10) yB! lOt 10) «F I! lOt 10) f PHI! 10,1 U) *GG( 10*10)* 
lOU! 10,10) ,0F I (10,10 ),F 120), WWI 150) 
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COMMON  £»C»btFI f PHI fGGtUQfQFI tF 

COMMON  Alt  lOflO)  f lOt  10)  tA3(  lOtlO)  t A<«(  10*10)  f Dl(  10»10)  » 

102(  10*10)  *Di(  10*10)  «0<t(  10*10)  *H(  20)  tCul  ( i0«10)«u02(  10*10) 

2«Al  1( 10*10) 

Common  w(2oO)«zz(1o*io)*li i lo) *ek( lo ) *o( io«io) *Q( lo* lo ) 

COMMON  0( 20*20) *TlMk*HS*HSTfcP*EPS»U 
COMMON  PFI  ( 10*  10)*0FH  10*101  *XINTA1(  10  * 10  ) *X  1NTA2  (10*10) 

1*YF 1 ( 10*10) ,OYF 1 ) 10*10) .YINr All  10*10) *YINTA2( 10*10) 

Common  k*m*mmk*  ieIo* inoex( lO) * jcounti lo) 

1»KCUUNT*IMETH*  NT*MTWU 
COMPCfcX^lo  Yl 10*10) 

CALL  CNVRT  l 1*Y  ) 

CALL  OIST 

IF  ( lMETH.bQ.2 )RfcTURN 
IF( iEiG.eo.o)on  To  10 
IM6TH=2 
RETURN 

10  00  17  I=1*K 

DO  17  J=1*K 
W( 1)=DREAL(C{ I* J)  ) 

W(2)=0IMAG(C( I* J)  ) 

IF  ( OABS  (»»'  ( n ) .LT.  1.0-1^)  Wl  I ) =0.00 
IF l DABS l Wl 2 ) ) .L  T.l.U-lA) W( 2 ) =0.00 
1/  C l I * J ) =OCMPLX  IW  l 1 ) , W(  <;  ) ) 

20  CALL  CNVRT  I2*Y) 

R£1 URN 

end 

SUBROUTINE  OIST 

IMPLICIT  RcAL*8l A-H) *R£AL*6(U-Z) 

CUMPLr:X«l  t>  £(10)*C(10*10)*B(10*10)*FI(10*10)*PH1(10*10)  *GG(  10*10)  * 
1QU(  10*10)  *CF  I (10*10  )*F(  20)*  nWdbO) 

COMMON  £*C*B*FI *PHI *GG*UO*QFI *F 

COMMON  Al( 10*10) *A2( 10* lu) *A3( 10*10) *AA( 10*10) *01( 10*10) * 

102(  10*10)  *03(  10*10)  *C^(  10*  10)  *H(20)  *OJl  ( 10*  10)*U04::(  10*10) 

2*A11( 10*10) 

COMMON  W(2o0)*ZZ(lU*i0)*El(lO)*EK(10)*G(10*10)*Q(10*10) 
common  o( 20*20) *time*hs*h step *EPS*U 
COMMON  PFI ( 10* 10)* OF  I ( 10  * 10 ) *X INTA I ( 10 * 10 ) * X INTA2 ( 10*10) 

1*YF1 ( 10*10) *0YF 1 ( 10*10) *Y1NTA1 ( 10*10) *YINTA21 10*10) 

COMMON  K*M*MMK.  I t I O* I NOE X ( 1 0 ) * JCuUNT ( 1 0 ) 

1»KCuUNT*IML IH*  nT*MTWO 
10=0 
101=200 
lOl ST=0 
00  20  I =1 *K 
JCOONT ( I ) =0 
FI( 1)=01MaG(E( I ) ) 

£R( I )=OKEaL(E( I ) ) 

C TEST  STABILITY  BY  CONSIDERING  ZEROS  OF  Al. 

IF( £R( I ) .LT.O.OOIGO  TO  20 
JcOUNT( I ) =1 

IF( £R( I ).£O.O.DO)GO  TO  20 
JCOONT ( 1) =2 
20  CONTINUE 

C TfcST  FOR  REPEATtO  ElGtNVALULS  OF  Al. 

IEIG=0 

IF(K.EO.l)GU  TO  3S 

KMlxK-1 

DU  30  J«1«KM1 
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JPl =Jt  I 
Du  30  I-JP1,K 
If  ((K(  J)  g.LkI  I ) ) IDs] 

U (t  1 ( J>  I ( I)  ) 1DI«  I 

U ( I 0. tu. I u( I 1 1 10= 11  10*1 
jO  CUNTINUt 
35  CONTlNUfc 
RfcTURN 
END 

SUBROUTINE  FUN') 

Implicit  re al  =^=8 i a-m), he al*b  10-2 1 

CUMPLEX«lO  E ( 10  ) tL(  10»IU»  tlM  10.  10  ) »FI  ( lOt  10)  f PHK  lOtlO  ) »G0(  10»  10  » » 
lUUUO.lO)  fUF  I ( 10.10  ) tF(20).HW(  150) 
cummon  e . c » b . f I . phi . go . 00 . of  1 . f 

CUMHON  A1(10.10).A2(10.10).A3(I0.10).A^(10.10)«01(10.10). 

1021  10.10)  .C3(  10 .10)  .D4(  10.10)  .H(  20)  .001  ( 10.10)  .twD^  (10.10) 

2. AIK  10.10) 

COMMON  W(200).2Z(10.10).EI(10).ER(10).G(10.10)«Q(10.10) 

CuMMON  0( 20.20) .TIML.HS.HSTEP.EPS.U 
COMMON  PFI ( 10. 10) .UF I ( 10. 10) .XINTAl ( 1 0 . 1 0 ) . X I NT A2 (10.10) 

1 .YF I ( 10,1 0) ,DYF I(in,10).YlNTAl(l0.lJ),YINTA2(10,lO) 

Cummon  k.m.mmk.  le i u. inofx( lo ) , jcount ( lo ) 

1 .KCUUNT . 1 M( 1 H.  NT.MTWU 

c create  a diagonal  matrix 

21 K0=0. 

OU  20  1=1, K 
FI(  l.I)=0.DU 

IF(  ER(  I ) .LT.0.D'J)G(J  TO  o 
IF ( E I ( I ) . LT .0.00)00  TO  6 

5 F M I . I ) =CUFXP)t  ( I X'T  IMF  ) 

CU  TO  7 

6 IF(COaBS(E( 1 X-TIME) .LT.179.U0)GU  TO  5 

7 IF (CUABS (E ( I ) ) .Eu.ZERO) GU  TO  10 
OFI ( I. I )= (FI ( 1. I )-l .00) /E( I ) 

GU  TO  15 

10  OFK  1.  1)  =TlMt 
15  00  20  J=l.K 

IF(J.EQ,I)GU  TO  20 
FI( I. J)=0, 

OFI ( I.J ) = 0. 

20  CONTINUE 

CALL  OCMPY( C.K.R.FI ,K .K.PHl ) 

CALL  OCMPY(PHI,K,R.B,K.K,bG) 

C COMPUTE  MATRICES  G AND  W 

CALL  OCMPY(C.K.K.OF I .KfR.PHI ) 

CALL  DCMPY( PHI.K.K.b.K.K.UO) 

IMG=0 

IMy=0 

00  50  J=1.K 
un  50  1=1. K 
0u=D1MAG( GG( I .J ) ) 

Ou=UIMAG(OU( I.J) ) 

IF(OABS(OG) .GT. 1.0-3)  IMG=IMG^1 
IF(0ABS(00) .uT.l.n-3)  IM0=IMQ*1 
50  CONTINUE 

1F( IMG.GT.O) WRITE (6.60) IMG 
1F( IM0.GT.0)WR1TE(6«  70) IMO 

60  FORMaT(3X.I2.*  of  THE  IMAGINARY  VALUES  OF  G MATRIX  HAS  A VALUE 
1 uREATER  THAN  lE-3' ) 
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70  F0KMAT(3Xf I2f • OF  IHt  IMAGINARY  VALUtS  OF  U MATRIX  HAS  A VALUE 
I GREATER  THAN  lE-3*  J 
OU  flO  J=ItK 
OU  80  I =1 *K 
G(  I f J)=UREAL(GG( I«J  t ) 

80  Q( I « J I =UReAL(QO( I t J ) ) 

RF  TURN 
ENU 

suhruutine  uMATKX 

IMPLICIT  KFAL*^! A-M) »KFAL«e(0-Z ) 

CuMPLt;<=>l6  t ( 10  I fC(  lOflul  *P  ( lOt  10  ) »FI  ( lOt  I0»  tPHK  10,10)  ,GGI  10,10) 
IQU( 10,10) ,QF I ( 10,10 ) ,F( 20), WW( 1^0) 

Common  E,c,b,Fi ,phi ,gg,co,ofi ,f 

COMMON  AI(I0,10),A<:(in,I0),A3II0,I0),AA(I0,10),DI(I0,10), 

1 02  < 10,10)  ,Oo(  10,10)  ,C<t(  IQ,10),H(?0),Q01(10,IO),U02(  10,10) 

2,A1 1( 10,10) 

COMMON  W(200),Z2(10,I0),EI(Io),Ek(I0),G(10,10),0(10,I0) 
common  0(20,20) ,TIML,HS,HSTeP,EPS,U 
COMMON  PFI  ( 10, iO),OFI ( IO,IO>,XINTAI( 10,10) ,X1NTA2( 10,10) 

1,YFI ( 10,10) ,OYFI ( 10,10) ,YINTA1( 10,10) ,Y1NTA2( 10,10 ) 

COMMON  K,M,MM)C,  I E I b,  INDEX  ( 10  ) , (COUNT  ( 1 0 ) 

1,KC0UNT  , IMC-TH,  NT,MTW0 
MP l = M*  I 
MP1  = M*  1 
MMK=M-K 
K1=K*1 

CALL  OMMPY( U,K,K,Ul ,K,K,U01 ) 

00  50  J=1,K 
00  15  1=1,K 
0( I,J)=G(I,J) 

15  01  1 *M, J) =QOl I I , J) 

00  2o  L=Kl,M 
I=L-K 

0(L,J)=TIME«(A3(I,J)*03(I,J)) 

0(L*M,J)=0. 

20  CONTINUE 

DO  30  1=1, MhK 

30  g01(J,I)=  A2(J,I)«02( J,l) 

50  CONTINUE 

CALL  DMMPYl  Q,K ,K,001 ,K,MMK,002 ) 

00  70  J=l,MMK 
00  bO  I=1,K 
L = J*K 

0( I,L)=  002( I,J) 

60  0(I*M,L)=0. 

00  70  I=K1,M 
L = J*K 
I I=I-K 

0( I ,L)  = TIME«OA(  II, J) 

0( I*M,L) =TIME#AA( I I, J) 

70  CONTINOE 

00  80  I=K1,M 
0( l,i)=0( I,I  )*I. 

80  CONTINUE 

00  100  I*1,MTM0 
00  90  J=MP1,MTM0 
90  D(I,J)=0. 

1F( I.CT.M)G0  to  100 
0( I, l*M)al. 
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100  LONTINUE 
RETURN 
END 

subroutine  MOUEIG 

IMPLICIT  REALOai A-H) #KEAL«fi(U-2 ) 

CuMPLEXOlb  k«  10  I tCUOtlO)  tBl  10,  10  J ,FI  ( lOf  10)  »PHI(  lOf  10>  fGG(  10,10)  t 
1Q0(  10,10)  ,OFI  ( 10,  10  ),EUO),RW(  150) 

COMMON  E,C,B,FI ,PHI ,GG,OQ,OFI ,F 

COMMON  Al(in,10),A2(10,10),A3(10,10),AA(10,10),01(10,l0), 

1021 10,10) ,03( 10,10) ,041 10,10) ,H( 20) ,0UI ( 10, 10 ),u02( 10,10) 

2, AIK  10,10) 

COMMON  WI200),ZZ(10,10),EI(10),ER(10),G(i0,I0),0(l0,10) 

COMMON  01 20,20) , T I Mfc ,hS ,HSTE P , F PS , U 
COMMON  PFK  10,  10),OFK  10,10  ) ,X1NTA1  ( 10, 10  ) ,X  INTA2  ( 10,10) 

ItYFU  lOtlO>tDYFl(  10,  lu  ),  y INT  A1  ( 10,10  ),yiNTA2(  10,10  ) 

COMMON  K,M,MMK,  IEIG, 1N0EX( 10) , JCOUNTI 10) 

1, KC0UNT,IMETH,  NT,MTWO 
C0MPLEX*I6  y( 10,10) 

M2=MTW0 
KC0UNT=0 
0NE=l.n0 
UO  10  1=1, M2 
10  1NUEX(1)=0 

CALL  CNVRT  n,y) 
no  20  1 = 1, M<; 

IFICOABSIFC 1 ) ).LE.ONE)GO  TO  20 
RCOuNT=KCOUNT»l 
INOtXIKcOUNT )=l 
20  CONTINUE 
RETURN 
ENO 

SUBROUTINE  LNVRT  II,y) 

IMPLICIT  RtAL<>rt(A-H),REAL<‘8tO-Z  ) 

CuMPLkX«l6  L(  10  ) ,C<  10,10)  ,P(  10,  10  ),F  I ( 10,10  ),PH1(  10,10)  ,Gu  I 10,10), 
1QU(  10,10)  ,UF  1 UU,  10  ) ,»-(20)  ,RW(  1 50  ) 
common  E,C,U,FI ,PHI ,GG,0Q,QFI ,F 

common  A1(  10,10),A2I 10,10), A3I 10,10 ),AA(10,10),D1( 10,10), 

1021  10,10) ,03(  10, 10) ,04 ( 10, 10) ,H( 20) ,0011 10, 10), 002 ( 10,10) 

2, A1 If  10,10) 

COMMON  W(200),ZZ( 10,10) ,EI ( 10) ,EK( 10) ,Gt 10,10) ,0( 10, 10) 

COMMON  Of  20,20)  ,TIME,HS,)<STeP,EPS,U 
COMMON  PFI t 10, 10), UFI f 10,10) ,X INT A If  10,10 ),X INT A2f 10,10) 
l,yFIf 10,10),DYFlf 10, I0l,yi NT Alf 10,10 ),YlNTA2t 10,10) 

COMMON  IC,M,MMK,  I E I O,  INOE  X f 10  ) , JCOUNT  f 1 0 ) 

1,KCUUNT,IMETH,  NT,MTWO 
COMPLEX^lb  YfK,K) 

IFf 1-2) 10,20,30 

10  call  EIGRFfAl,  K,10,2,E,C,10,W,IER) 

WRlTEfb,  15)  ItP. 

15  format  1 2 X,* ERROR  PARAMETER  FROM  EIGENVALUE  ROUTINE  HAS  THE  • 

1, 'VALUE  *,13,/) 

IFf Wf 1) .GE.1.00)G0  TO  17 
WKITEf6,16) 

16  FORMATf 2X,»THE  tIGFNVALUE  ROUTINE  PcRfOkMeO  MELL*,/) 

GO  TO  40 

17  IFf Mf 1) .GT«100.00)GU  TO  19 
«fRITEf6,18) 

Id  FORMATf 2X,* the  EIGENVALUE  ROUTINE  PERFORMED  SATISFACTORILY*,/) 

GO  TO  40 
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1 

19 

WKI rEI6t21  ) 

t 

21 

FORMAK 2X, • THE  EIGENVALUE  ROUTINE  PERFORMED  POOKLY't/) 

0 1 

GU  TO  40 

20 

00  2 IIsl,K 

1 

00  2 J = UK 

i 

2 

Y( II«J)=C( iltJ) 

u-1. 

\ 

CALL  MA23B0(Y«KtKf WW«U  ) 

’ 

IFtU.OE.O. 00)00  TO  3 

IMETHs? 

i 

WRI TE (6y8U) 

80 

FORMAT) 2X, ‘MATRIX  INVERSION  ROUTINE  HAS  PRODUCED  UNRELIABLE  RtSuLI 

lS»./«2Xy'THEREFURh  NUMFRICAL  INTEGRATION  METHOD  WILL  dE  USEO.'t/) 

j 

\ 

GO  TO  40 

3 

DO  5 II=ltK 

OU  5 J=itK 

5 

B( I I* J)=Y( I I ,J) 

GO  TU  40 

30 

Call  EIGRF  ir)»MTWUf20»0yFf2Z»I0yWtIER) 

WRI rE(6yl5) lER 

40 

RETURN 

END 

» 

SUURUUTINE  DMMPYCDPMl ,MMM*NNNI t DPM2 t NNN2 t LLL  t DPRES ) 

c 

matrix  multiplication 

» 

c 

DPMI  — input  is  an  MMM  X NNNi  MATRIX 

• 

c 

PUT  Is  UROEREO  BY  THE  COMPUTER  AS  A 10«10 

f 

c 

0PM2  — input  IS  A NNN2  X LLL  MATRIa 

1 

c 

BUT  IS  uRUEREU  OY  THE  COMPUTER  AS  A 10«10 

r 

c 

OPRES  — OUTPUT  IS  AN  10  X 10  SOUaRE  MATRIX 

1 

c 

iM. 

n.  INPUT  MATRICES  UF  uRuER  A 

i 

c 

IMPLICIT  REAL=>8(A-H)  ,REAL#fl»0-Z  ) 

tf 

I 

1 

r 

DIMENSION  DPMI)  lOylUl fOpN2( lOy ID > yOPRcS ( lUtlO) 

1 

c 

CHcCK  THAT  OlMtNSIONS  UK 

1 

IF(NNN1.NE.NNN2I  GO  TO  70 

• 

f 

c 

1 

OU  SO  K-I*LLL 

j 

c 

OUTSIDE  LOOP 

1 

2 

DO  50  1=1, MMM 

f 

c 

INSIOt  LOOP 

y 

c 

1 

c 

BUILD  EACH  ELEMENT  OF  THE  PRODUCT 

1 

CTEMP=0. 

■ 

00  40  J=lyNNNl 

c 

COLUHN«RUW 

CTEHP=CTEMP  ♦ DPM H 1 , J ) ^OPMZ ( J, K ) 

40 

CONTINUE 

• 

c 

r 

50 

OPRES) I ,K)=CTEMP 

• 

c 

RETURN 

c 

INPUTS  NUT  EUUAL 

V 

70 

WRITE(6t201 )NNNi,NNN2 

201 

FORMAT) IX, • UMMPY  DIMENSIONS  FED  lU  NOT  COMPATIBLE  FOR  *, 

1 /,*  multiplication  — EXECUTION  TERMI MATED* ,/ ,2X ,• NNNI *• , 

■ t 

i I 1,5X, *NNN2=*, 13) 

no  nno  o o oo  nor>r»onno 
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STuP 
bNO 

5UUR01JTI  Nt  Ot.MPY(f.UMPLl  tMM.MtNNM  t CUMPL^  ♦ NNi<2  . LLL  tC  MPKP  :» » 

COMPLEX  MATkIX  MULTIPLILATIJN 

COMPLl  — INPUT  15  AN  MMM  X NnNI  LOMPlFX  MATRIX 

BUT  IS  UPUEKEU  BY  THE  COMPUTER  AS  A 10«10 
CUMPL2  — INPUT  IS  A UNNZ  X LLL  COMPLEX  MATRIX 

BUT  IS  ORUEREU  oY  THE  COMPUTER  AS  A 10*10 
CMPRES  — OUTPUT  IS  AN  10  X 10  SQUARE  COMPLEX  MaTKIx 

N.  B.  INPUT  MATRICES  OF  ORDER  I 

COMPLEX#  16  COMPLK  10.  lO  ) .COMPLZ  ( 1 0. 1 0 ). CMPRES  11  0.  1 0 ).  CT  EMP.LMPLX 

CHECK  that  dimensions  OK 

IFINNNI.NE.NNNZI  GU  TO  70 

UO  SO  K=i.LLL 
OUTSIDE  LUOP 
DO  SO  1=1. MMM 
INSIDE  LOUP 

BUILD  EACH  ELEMENT  UF  THE  PKOOUCT 
CTEHP=CMPLX(0..0.) 

DU  <tO  J=l.NNNl 

COLUMNOROW 

CTEMP=CTEMP  ♦ CUMPLK  I.  J»#COMPLZ(  J.KI 
AO  CONTINUE 

SO  CMPKESl I.K) =CTEMP 
C 

RETURN 

C 

C INPUTS  NOT  EQUAL 

70  MRITE(6.201 INNN1.NNN2 

201  fokmatiix.*  ocmpy  dimensions  fed  In  NOT  compatible  for  ». 

1 /.'  MULTIPLICATION  — EXECUTION  T EKM I NATEO • ./ . 2 X . • NNNI = • , 

2 I3.SX. ‘NNNZs*. 13) 

STUP 

ENO 

SUBROUTINE  NUMERIC IPASS. IN. HINT ) 

IMPLICI T REAL#8( A-M) .KEAL#8C0-Z ) 

COMPLEX#! 6 E(  1C ) .C< 10.10) .D( 10.10 ).F I ( 10.10) .PHI ( 10.101 .GG( lO.lO) . 
1QU(  10.10).OFI(10.10).F(20).kWnSO) 

COMMON  E.C.b.FI .PHI .GO.QQ.OFI.F 

COMMON  Al( 10.10). A2< 10 .10). A3 ( 10.10 )«AA( 10.10). OK  10.10). 
lD2t 10.10) .03 ( 10.10) .0A( 10.10 ).H( 20) .QOK 10.10 ).Q02( 10.10) 

2.A1 1( 10.10) 

common  W(200).ZZ(10.10) .ET ( lOl.ERI 10).b( 10.10) .Qt 10.10 ) 
common  0«  20.20) . TIME. HS.HSTEP. EPS. U 
COMMON  PFldOf  LOlfUFlt  10.  10  ) .X  I NT  A 1 f 1 0. 10  ) » X 1 NTA2  ( 10*10) 
l.YFll  10.10).OYFII10*10)*Y1NTA1(  10  * 10  ) * Y INT  A2  ( 10*  lo) 

COMMON  K.M.MMK*  I E 1 U. INDEX ( 10 ) * JCUUnT ( 1 0 ) 

1.KCUUN7.1METH*  iiT.MTWU 
C PERFORM  NUMERICAL  INTEGRATION  PROCESSES  To 
C DETERMINE  MATRICES  G AND  U. 

DIMENSION  UNIT (10* 10) 

IF( IPA$$.NE.O)GO  TO  40 
C initialize  VARIABLES  FUR  SINT. 
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C calculate  STEP-LtNGTH  FOR  NUMERICAL  INTEGRATION  INTERVAL. 

LSTEPsO 
TAU^O.OO 
OU  10  I=l,K 
00  5 J=1«K 

5 UNI T( It J) -0.00 
10  UNIT! If  11=1.00 
Ou  15  I=l,K 
OU  15  J=l.K 
PFnitJ»=UNIT(I,J» 

XlNTAUlf  JI=UNir(ltJi 
X1NTA2I It J)=0.00 
0FI(ItJ)=Ali( It Jl 
YFI ( I tJ 1=0.00 
OYFI  ( 1 1 J) =UNIT( ItJ) 

YINTAU  It  J)=0.00 
15  YINTA2( It J»=0.00 

40  IFILSTEP.GT.OIGO  TU  45 
CALL  STEP  I INtHlNTtLSTEP) 

45  NAM=-l 
IPASS=l 

tau=tau*hs 

50  IPASS=-1PASS 

CALL  SINT ( IPASStNAM) 

IFI NAM) 52t52t54 

52  Call  OMMPYIAUtKtRtPFltKtRtOFl) 

GO  TO  60 
54  OU  56  I=ltK 
00  56  J=ltK 
56  OYF  I ( It J) =PFI ( 1 1 J) 

60  IFI IPASS.lt. 0)G0  TO  50 
NAM=-NAM 

IF (NAM.uT.OIGO  to  50 

tau=tau»hs 

IF(TAU.LE.TIME)GO  TU  50 
TAU=TAU-HS 
DO  70  I=ltK 
00  70  J=ltK 
Gl  I tJ)  = PF  I( 1 1 Jl 
70  01  ItJ)=YFH  I.J) 

RETURN 

FNO 

SUoROOriNE  STEPI ItHlNTtLSTeP) 
implicit  RCAL*8(A-HltREALP8(0-Z I 

C0MPLEX016  El  LO)tCI  lOtlO) t8( lOt 101 tFIl lOtlOltPHK iOtlU) tGG( lOt lul t 
lUUl lOtlO) t0FI(I0tl0)tFI20)tWMII50) 

Cummon  e t c 1 1>  t f 1 1 phi  t gu 1 00 1 of  1 1 f 

COMMON  All 10tI0ltA2l lOtlOltAil 10tlOltA4<IOtlO)t01(10tlO)t 
1021 lOtlO) t0j{ 10,10) t04l 10, 10 ) fH I 20 ), 001 1 10, 10) ,002  I 10 1 10) 

2, Alll 10,101 

COMMON  WI2uO),ZZllOtlO)tEI I lU ) ,ERI 10 ) , Gl 10, 1 0 ) ,0 1 1 0 1 10 ) 

COMMON  0120,20) ,TIMEtHStHSTEP, EPS, 0 
COMMON  PFll 10tl0),0Fll 10,10) ,X1NTA1I 10,l0)tXlNTA2l 10,101 
1, YFI I 10,10) tOYFl 1 10, lu), YINTAU 10,10) ,Y1NTA2 1 10,10) 

COMMON  K,MtMMKt  Id  Gt  INuEX I 10 ) , OCOUNT  I 1 0 ) 

1,KC0UNT,IMETH,  NT,MTWO 

C establish  step  LENGTH  FOR  NOMcRlCAL  INTEGRATION  FuR 
C STABILITY  ANO  ACCURACY  OF  NUMERICAL  SOLUTION. 

C SUPPLY  STEPLENGTH,HStFOR  EACH  INTERVAL  I HI  I ♦! ) tHI I ) ) . 
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1F(  UGT.ZIGJ  TO  bO 
IK  I.GT.I  »GU  in  20 
Hri=H(  1 1 
GO  TO  21 

20  HM  = H(  n-H(  I-  U 
21  XVAL=  HH  /HSTcP 

NSTEP  =0If<T<XVAL)  ♦! 

HS  = HH  /uFtOAHNiTfcP) 

GO  TO  60 

50  in  HiNT.Lfc.O.  )G0  TO  20 
LSTLP=1 
60  KfcTURN 
END 

SUbRQUTI.NE  SINT  ( IPASitNAM) 

IMPLICIT  RCAL*8(A-H»fKEAL*8IO-Z I 

C0MPLEX«I6  E I 10)  «C(  iOf  10)  yCI  I0«  I0)*FI  ( 10*10)  tPHK  10*10)  tGo(  10«Iul 
1QU(  10*10)  *QF  I ( 10*10  )*F(  <^0)  * MW ( 150) 

COMMON  E * C * b * F I * PHI  *GG*  * OF  I * F 

COMMON  Al(  10*10)  *A2(  10*  10)  *Ai(  10*10  )*A<i(  10*10)  *01  ( 10*10)  « 

1D2( 10*10) *03 ( 10* 10) *0^( 10*10 )*H(20) *001 ( 10* 10)*u02l 10*10) 

2*A11( 10*10) 

COMMON  W ( 200  )*ZZ(  10*10)  *EI  ( 1 0 ) * EK  ( 1 0 ) * ( 1 0*  1 0 ) *0  I 1 0*  10  ) 

COMMON  D«  20*20) * T IMc *HS *HST EP* E PS *U 
COMMON  X IN [( 10*10) »DXI NT ( 10*10 ) *XINTA1 1 10* 10 ) * XINTAZ 1 10* 10 ) 
l*yiNT( 10* 10) *DY1NT( 10 * 1 0 ) * Y 1 NT A 1 ( 10  * 10 ) * Y 1 NT A2 ( 10*10) 

COMMON  K*M*MMK*  IfclG* INDEX! 10 ) * JCOUNTI 10) 

I*rCUUNT*IMFTH*  NT*MTWU 
HaN2=HS/2.00 
IF(IPASS)1*I*3 
1 IF(NAM.GT.O)GO  TO  15 
on  4 I =1 *K 
DO  4 J=1*K 

XINT(I*JI  = XINTAI(I *J)>HS«DXINT ( I * J) 

4 X1NTA2(  I «J)=:DX1NT(  1*J) 

GU  10  10 

15  DO  16  I=1*K 

DO  16  J=1*X 

YINT ( I ,J )=YINTA1( I *J) ♦HSOOYINT ( 1 *J) 

16  Y1NTA2(  l*J)-nYINT( I *J) 

GO  TO  10 

3 IFINAM.GT.OIGO  TO  25 

DU  7 I=1*K 
UO  7 J=1*K 

XINTI I «J)3XINTA1( I * J)»HQN24(XINTA2( I* J) ♦OX INTI 1*J  ) ) 

/ XI.^TAI  ( I *J)sXlNT(  I *J) 

30  TO  10 
25  00  18  I=1*K 

DO  Id  J=1*K 

YINT I I*J)=YINTA1II*J)*HON2«IYINTa2I I*J)«dYINT<1*J)) 

18  YINTAII 1*J)=Y1NT( 1*J) 

10  CONTINUE 
RETURN 
ENU 

SUBROUTINE  LAMBDA 
IMPLICIT  REAL«8IA-H )*REaL*8(0-Z ) 

C0MPLEX«16  E( 10) *C< 10«lu) *H( 10*10 )*FI( 10*10)*PH1( 10*10) *GG( 10*1J) 
lOUl 10* 10) tOFI 110*10 )«F(20) *MW( 150) 

COMMON  E*C*b*FI*PHI*GG*wO*OFI*F 

COMMON  Ail  10*10 )*A2( 10* 10) *A3I  10*10 )*A41 10*10) *D1 1 10*10) * 
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1021  lOt 1 01 «03( iOt 1 j) tUHl i0«l0 t «H(2UI f QUl ( lUt iO ) *w02 ( iC « 1 J ) 
2«A11(  iOtlOl 

CUMHON  W(?00)  tZZ(  lUt  10)»t:i(10)«ER(i0)«o(l0tl0)t0(l0ti0> 
COMHON  0( 20,20) , T t , HS* HST fcP , E H S ,U 
COMMON  pm  lOf  iO),uFi  ( lOf  lOfXiNTAK  1 0 , 1 0 ) , X IN  TAii ) 10 , 1 0 ) 
UYFK  10,10)  ,DYFI  ( 1C  , 1 0 ) , Y 1 NT  A 1 ( 10  , lo  ) » Y 1 NT  A2  ( 10 , 1 0 ) 
COMMON  K,M,MMK,  I t I O, 1 NOEX ( iO ) , jCuUNT( 10 ) 

1,KC0UNT,IMETH,  NT,MTWO 
DIMENSION  HNO(IO) ,HnI( 10) 

NTi=NT-l 
HOSrK=H( I ) 

00  10  1=1, NT! 

HSTAb=K< I »1 )-H( I ) 

IF IHSTaB.LT.HUSER) HUSER=HSTAd 

10  continue 

WRITEI6,15)HUSEK 

15  F0PMAT«2X, 'MINIMUM  USER  STEPLtNCTH  INTERVAL  IS  ',1:13. t>) 

£PS=0. 100 

IF( ER( 1 ) .cU.O,DO)cn  TO  16 
HMIN  =(-2»U0=>EPS)/tRll  ) 

GO  TO  18 

16  HMIN=2000.00 
lb  DO  100  I=1,K 

IF(EK)  1)  .Ew.O.DOX'.U  TO  20 
HNO( 1)=I -2.D0«£PS) /Ek( I ) 

CO  Tc  30 

20  HN0<I)=10.c0 

30  IFCEK  I)  )100,A0,50 
<,0  HNI(I)=I0.0U 
CO  Tu  70 

50  IKHNO(  I)  .Ey.  10.Cj)c0  TO  oO 
OMB=  l-a.OORLR  ( I )<=riNU«  I ) ) «=<=0,25 
VH=»E1(  iX'HNOd  ) )/rMB 
IF(  VH.GE.0.<:00)c0  Tu  oO 
HLAM=HNO( I ) 

CO  TO  90 

60  OM=-kR ( r )/( E I<  I )#*H ) 

PSI=d.OO<=l  I2«00<'LPS)<=<''t) 
cXPN=l .OU/J.OO 
HNK  I ) =(OM«PSI  )R-r:XPN 
70  HLAM=0M1N1 ( HNO( 1 ) ,HnI  ( I ) ) 

90  1F<HLAM.LT.HMIN)HM1N=HLAM 
100  LUnTl.'JUE 

MSTEP=OMlNllHMiN,HUStR) 

WRITE! o, 105 IHSTEP 

1U5  FORMAT (2X,'N. I . STEPUENGTH  FROM  S/K  LAMBDA  IS  ',£13. 6) 
RETURN 

tNO  / 

subroutine  OUTPUT  1XX,L,N,X) 
implicit  RLAL<‘3(A-H),REAL«8(0-Z  ) 
dimension  DcMMYI20) ,XXI10,1U),XIL,N) 

Do  2 1=1, L 
00  2 J=1,N 
2 XII,J)=XXII,J) 

00  10  1=1, L 

WR1TE(6,5) (X( I , Jl, J= I,N) 

5 FORMATdH  , 10(F9.2,  2X  ) ) 

10  CONTINUE 

WHITE! 6,15) 

15  format  ! I'a,  lOO!  IH.  ) ) 

RETURN 

END 
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APPENDIX  IV 

LIST  OF  VARIABLES  USED  IN  SOFTWARE  PACKAGE 


Al, 

Dl 

A2, 

D2 

A3, 

D3 

A4, 

D4 

All 

M 

K 

KK 

KMK 

MKMK 

KTTWO 

H 

NT 

EL 

HINT 

IMETH 


E 

ER 

El 

C 

B 

W 

JCOUNT 


This  family  of  real  matrices  form  the  analogue  and  digital  pairs 
which  describe  the  partitioning  scheme  being  tested.  The  row 
and  column  dimensions  in  the  calling  program  are  10  for  each 
matrix.  The  actual  size  of  each  matrix  is  defined  by  the  user 
for  his  problem. 

Stores  A1  which  could  be  lost  during  execution  of  the  IMSL 
eigenvalue  routine  EIGRF. 

Number  of  first  order  d.e.'s  in  the  system,  i = Aj^,  once  reduced. 
(M  > K+1). 

Actual  row  or  column  size  of  Al,  D1  (K  10). 

= K X K;  number  of  elements  in  Al,  Dl. 

= K x (M-K)  number  of  elements  in  A2,  D2,  A3,  D3  (M-K  ^ 10). 

= (M-K)  X (M-K)  number  of  elements  in  A4,  D4. 

2M  (2M  < 20) . 

Real  vector  of  length  NT  containing  the  user  values  of  step- 
length  for  testing  the  partitioning  scheme. 

Length  of  the  vector  H. 

First  element  of  H;  defined  when  the  required  step- lengths  are 
at  equal  time  intervals. 

Time  interval  between  successive  step- lengths  H.  and  when 

EL  is  defined. 

= 1 The  eigenvalue  method  is  used  for  computing,  4>,  G,  Q. 

= 2 The  numerical  integration  method  is  used  for  ♦,  G,  Q. 

The  user  defined  value  for  IMETH  can  be  overwritten  by  the 
program  depending  on  the  value  assigned  to  lEIG. 

Complex  vector  of  length  K containing  the  eigenvalues  of  Al. 

Real  vector  defining  the  real  part  of  E. 

Real  vector  defining  imaginary  part  of  E. 

Complex  matrix  of  order  K x K of  eigenvectors  of  Al. 

Inverse  matrix  of  C. 

Working  space,  dimensioned  at  least  (2  > K)  x K,  used  in  the 
library  routine  EIGRF. 

Vector,  of  length  K,  which  is  used  to  print  the  results  of  the 
stability  analysis  of  G. 
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1ST 

IMST 

INST 

IBIG 


TIME 

FI 

QFI 

PHI 

GG 

QQ 

SOTE* 


TAU 

EPS 

HNO 
HNI 

HUSER 

IMIN 


Counters  for  recording  whether  each  eigenvalue  of  A1  has 
negative,  zero  or  positive  real  part.  The  result  defines  the 
stability  of  G. 

Counter  for  accumulating  the  number  of  repeated  eigenvalues 

lEIG  = 0 means  the  X.  are  distinct. 

1 

Set  to  for  i * 1,...,NT. 

Complex  matrix  defined  by  (exp(X ^h) , . . . ,exp(Xj^h))  . 

. ,exp 


(Xjh)-l 

Complex  matrix  defined  by  [exp — 


1 


NAM 

= 

-1 

s 

♦ 1 

I PASS 

- 

0 

3 

-1 

S 

♦ 1 

Complex  matrix  defining  (i)  CxFI  when  computing  G(h)  and 

(ii)  CxQFI  when  computing  Q(h). 

PHI  X B for  (i)  above. 

PHI  X B for  (ii)  above. 

GG  and  QQ  are  complex  matrices.  It  is  expected  the  imaginary 
part  of  these  matrices  will  be  less  than  10”*. 

FI ,0FI ,PHI ,GG  and  QQ  are  all  defined  in  subroutine  FUND  for  a 
range  of  step-lengths  defined  by  H.  The  order  of  each  matrix 
is  K,  while  the  row  and  column  dimension  in  the  calling  program 
is  10. 

Real  matrix  of  order  K;  = GG  for  IMETH  = 1 

= <b(h)  for  IMETH  = 2 

Real  matrix  of  order  K;  = QQ  for  IMETH  * 1 

.h 

<l>(T)dr  for  IMETH  = 2 
Both  G and  Q are  defined  for  a range  of  step- lengths  H. 


Set  in  Main  for  initialisation  in  integration  routine  SINT. 
Set  in  subroutine  NUMERI  for  operation  of  the  "predict" 
step  of  SINT. 

Set  in  S/R  NUMERI  for  operation  of  the  "correct"  step  in 
SINT. 


Keeps  tally  of  time  from  0 to  H^  for  numerical  integration. 

Describes  the  accuracy  of  the  numerical  solution  for  ♦(h) . 

Real  vectors  of  length  K,  defining  values  of  N.I  step-length  as 
functions  of  X^,  for  j «>  1,  ....  K. 

Minimum  step-length  interval  between  elements  of  the  vector  H. 
Minimum  of  HNO(j),  HNI(j)  over  the  j eigenvalues,  X^ . 


I 


i 
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HSn;l‘ 


US 

NSTEP 

XINT 


YINT 
DXINT  ( 

dyintJ 

XINTA2'; 
Y1NTA2 j 

XINTAl  I 
YlNTAl J 

PFI 

YFI 

DFI 

DYFI 

ZZ 

QDl 

QD2 

D 

F 


Mlnlmiini  of  IIUSliR,  IIMIN,  und  defined  ns  the  niimerical  integration 
step- length  In  subroutine  I^DDA. 

N.I  step- length  based  on  HSTEP,  but  modified  to  suit  the  interval 
(H.,  H.^j)  for  i = 1 (NT-1). 

Integral  number  of  steps  of  length  HS  in  the  interval 

(Hi. 

Variable  computed  in  numerical  integration  routine  SINT, 

where 

x(tn^l)  = ^ "Predict"  step 

and 

x(tn^j)  *(^n^  * h/2(x(t^)  + "correct"  step. 

Used  for  computing  G(h). 

Analagous  to  XINT,  but  used  for  computing  Q(h). 

Define  x(tj^)  in  SINT,  for  G(h),  Q(h)  respectively. 

Hold  previously  computed  values  for  DXINT,  and  DYINT 
respectively  from  the  "predict"  step  in  SINT. 

Hold  previously  computed  values  for  XINT  and  YINT 
respectively  from  the  "correct"  step  of  SINT. 

Defined  in  NUMERI  as  and  is  equivalent  to  XINT  from 

SINT, 

Defined  in  NUMERI  as  and  is  equivalent  to  YINT  from 

SINT. 

Defined  in  NIMERI  as  ‘**(fj^)  = Al.  ***(^n^ ' equivalent  to 

DXINT. 

Defined  in  NUMERI  as  Q(fjj)  = G(t^),  and  is  equivalent  to  DYINT. 

K X K working  matrix 

Real  K X K matrix  defined  by  Q x Dl. 

Real  K x (M-K)  matrix  defined  by  Q x (A2  D2) . 

Real  2M  x 2M  matrix  which  is  a function  of  matrices  G,  Q,  QDl, 
QD2,  A3,  D3,  A4,  D4.  Defined  in  Section  3.2. 

Complex  vector  defining  2M  eigenvalues,  of  matrix  0. 
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KCOUNT 

INDEX 


Counter  defining  the  number  of  with  modulus  greater 
than  unity. 

Vector  for  storing  the  value  of  those  with  modulus 
greater  than  unity. 
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APPENDIX  V 

SOLUTION  TO  THE  EQUATION  i(t)  = 6^(1)  + Cu(t) 

Consider  the  general  equation 

i(t)  = B;c(t)  + Cii(t)  (V.l) 

and  let  X{t)  be  a fundamental  matrix  of  the  n x n system 

i(t)  = B^Ct), 

i.e.  X is  a matrix  such  that  its  n columns  consist  of  n linearly  independent 
solutions  to  the  system.  Then,  by  definition,  X(t)  satisfies  the  matrix 
differential  equation 


ciX(t) 

dt 


BX(t) 


CV.2) 


A fundamental  matrix  is  always  non-singular  (by  definition),  and  is  uniquely 
determined  for  a given  set  of  initial  conditions (ref , 4) . 

We  must  prove  that  the  solution  for  (V.l)  is 


Jl(t) 


X(t)  X'' (0) 


X(t)  X’'(0  Cii  a) 


with 


* X(t)  X*'  (0)  fy, 


complementary 

solution 


l£,(t)  , say 


particular 

solution 


)C(0)  = fi. 


To  establish  this,  we  need  only  prove  that  the  particular  solution,  satisfies 

equation  (V.l),  since  the  complementary  solution  is  well  established(ref.4) . 

That  is,  we  have  to  show 


= B )c(t)  + Cii(t) 


(V.3) 


Put 


v(t) 


X(t)  Z(t) 
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Then 


i = xz  + xz 

= BXZ  + XZ  from  equation  (V.2) 

= Bii  * XZ 


We  have  from  equation  (V.3) 


we  must  have 


i = Bx.  + Cii 


XZ  = Cji 


1 .e. 


Z = x"'  (t)  Cii(t) 


z(t)  = r x‘' (o  cjiicn 

^ o 


and 


x:(t)  = f X(t)  X"' {{)  CjArt)  d{ 

j o 

Choose  X(t)  such  that  X(0)  =1.  By  definition,  since  X(t)  is  unique, 

Bt 


Also,  by  definition. 


and  hence 


X(t)  = e^ 


X‘'(t)  = e'*^  = X(-n 


X(t)  X‘*  (I)  - » X(t  - f) 


Thus 


t 

■ / 

J o 


^(t)  . / x(t  - n Cttcn  dt 

' o 
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and 


4(t) 


X(t)  *,  ♦ / X(t  - t)  Cii(t)  di 


If  ll(5)  lijij  ® constant,  say,  and  t = h 

h 

X(h)  = X{h)  io  * r x(h  - n d{  Cju^ 

o 


Replacing  h by  r,  it  can  be  shown  that 


i(h)  = X(h)  ip  + ! X(T)  dr 
" o 

However,  any  fundamental  matrix,  ‘I’(t)  , for  tha  system,  can  be  written  as  the 
product  of  the  unique  matrix  X(t)  and  a non-singular  constant  matrix  C. 


i.e.  d>(t)  = X(t)C 


Thus 


C = '!«■'  (0) 


and 

X(t)  = ‘^(t)  4>''  (0) 


giving 

.h 

i(h)  = ‘h(h)  <J>’'(0)  Zfl  + I 'Hr)  ‘•’■*10)  dr 


APPENDIX  VI 

RESULTS  OF  APPLICATION  5.1 


IIC««VAi.Uf5  01  >1  (C0IIPL1.X  eLSIlEllS) 


■0«  « 0.0  0.J«02O-01  -0.64270-02  0.99290*00  0.203S0»00  0.0  0.0  0.0  0.0 

0.10000*01  0.0  0.0  0.0 
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2.  Flowchart  of  the  software  package 


WRE-TR-1882 (A) 
Figure  5 


Figure  5.  Block  diagram  of  height  control  system  simulation 
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Figure  6.  Partitioning  scheme  for  height  control  system  simulation 
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Figure  7,  Partioning  scheme  for  tracking  radar  simulation 
with  instantaneous  analogue  response 
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Figure  8.  Stability  boundary  for  the  p “edictor'cerrector 
numerical  integration  method 


